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Linear  feedback  controllers  and  estimators  have  been  designed  from  the  governing  equations  of 
a  channel  flow,  linearized  about  the  laminar  mean  flow,  and  a  layer  of  heated  fluid,  linearized  about 
the  no-motion  state.  Spectral  decomposition  involving  a  two-dimensional  Fourier  expansion  and 
a  Chebyshev-Galerkin  projection  cast  these  linearized  equations  into  state-space  form  that  decou¬ 
pled  to  independent  Fourier  wavenumber  sub-systems.  The  control  law  are  designed  by  applying 
Linear  Quadratic  Gaussian  (LQG)  synthesis  to  these  sub-systems.  The  size  of  the  controller  is 
reduced  by  both  limiting  the  number  of  sub-systems,  on  which  LQG  synthesis  is  applied,  as  well 
as  applying  system  theoretic  model  reduction  techniques  to  each  sub-system.  This  methodology 
has  produced  highly  effective  controllers  to  suppress  convection  in  a  heated  fluid  layer,  but  has 
found  only  moderate  success  with  the  channel  flow.  While  the  Oberbeck-Boussinesq  equations 
(heated  fluid  layer)  provides  a  direct  measure  of  Rayleigh-Benard  convection,  the  Poiseuille  flow 
equations  do  not.  The  feedback  control  laws  for  channel  flow  could  only  indirectly  affect  viscous 
drag.  An  open- loop  optimization  has  been  applied  to  the  channel  flow  control  problem  in  an  effort 
to  capture  more  of  the  nonlinear  dynamics  and,  thereby,  affect  the  viscous  drag  directly.  During 
these  experiments,  it  has  been  discovered  that  upstream  traveling  waves  of  blowing  and  suction  not 
only  reduces  the  skin-friction  drag  in  a  channel  but  also  sustains  it  below  laminar  levels. 
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1  Introduction 


The  potential  benefits  of  controlling  flows  to  reduce  drag  are  significant.  Commercial  airlines 
consume  about  1.5  billion  barrels  a  year.  A  flow  control  scheme  that  reduces  drag  by  20%  could 
save  $15  billion  per  year  (at  $50  per  barrel).  Reduced  drag  can  also  provide  faster  speed  instead 
of  savings  in  fuel  consumption.  In  the  past  decade,  flow  control  has  shifted  from  depending  on  the 
investigator’s  intuition  built  on  empirical  observations  to  more  systematic  approaches  focusing  on 
the  dynamics  of  the  flow  field  [1,2].  In  our  studies,  we  have  applied  modern  control  theories  to 
linearized  Navier-Stokes  equations  with  suction/blowing  actuation  along  the  channel  walls  with  no 
net  mass  transfer.  The  resulting  controllers  were  implemented  and  evaluated  in  a  direct  numerical 
simulation  (DNS)  of  a  turbulent  channel  flow. 

Spectrally  decomposing  the  Navier-Stokes  equations  linearized  around  the  laminar  profile  in 
a  channel  results  in  a  dynamical  model  in  state-space  form  that  decouples  to  independent  sub¬ 
systems  by  wavenumber.  A  parallel  architecture  of  compensators  in  wave  space  is  thus  possible  [3] . 
Control  laws  are  synthesized  for  only  certain  wavenumbers,  reducing  the  size  of  the  overall  con¬ 
troller.  Furthermore,  the  individual  sub-systems  can  be  reduced  using  modern  control  theories 
(based  on  observability  and  controllability).  Although  the  linearized  equations  are  missing  the 
nonlinear  terms  in  the  Navier-Stokes  equations,  Linear  Quadratic  Gaussian  (LQG)  controllers  of 
just  a  few  wavenumbers  have  been  moderately  successful  in  reducing  the  viscous  drag.  We  have 
reported  approximately  15%-18%  reductions  in  previous  reports  [4].  However,  further  reductions 
have  been  difficult  to  achieve. 

The  primary  obstacle  to  improved  performance  appeared  to  be  the  lack  of  a  cost  function  in  the 
linearized  equations  that  directly  relates  to  mean  skin-friction  drag.  The  same  control  methodol¬ 
ogy  applied  to  the  suppression  of  Rayleigh-Benard  convection  in  a  heated  layer  of  fluid  produced 
remarkable  results  [5,  6,  7,  8]  (Appendix  A,  B,  and  C).  The  Oberbeck-Boussinesq  equations  pro¬ 
vided  a  direct  quadratic  cost  of  convection;  however,  the  flow  control  problem  only  had  indirect 
cost  functions.  This  motivated  an  investigation  into  periodic  optimization  algorithms.  By  inducing 
small  perturbations  of  the  open-loop  controls  and  periodic  flow  field  conditions  in  the  nonlinear 
DNS,  direct  numerical  gradients  associated  with  the  drag  were  computed.  These  gradients  were 
used  to  make  improvements  to  the  control  history  and  initial  flow  conditions  to  directly  improve 
drag  reduction. 

In  the  process  of  simplifying  this  procedure,  we  discovered  that  wall-normal  blowing  and  suc¬ 
tion  on  the  upper  and  lower  walls  in  the  form  of  two-dimensional  waves  traveling  upstream  could 
reduce  the  skin -friction  drag  to  below  laminar  levels.  Defined  as  simple  streamwise  sine  waves, 
these  control  waves  do  not  change  the  net  mass  flux  of  the  channel.  Even  more  remarkably,  this 
behavior  appeared  to  be  a  sustainable,  steady-state  solution.  Although  the  properties  of  the  nonlin¬ 
ear  dynamics  that  results  in  the  sub-laminar  drag  are  not  completely  understood,  we  show  that  the 
Orr-Sommerfeld  equation  is  a  good  predictor  of  the  nonlinear  channel  flow’s  response  to  traveling 
waves. 
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2  Achievements 


The  linearized  three-dimensional  Navier-Stokes  equations  as  the  basis  for  designing  fully  three- 
dimensional  reduced-order  controllers.  Improvements  to  the  estimation  in  the  LQG  synthesis 
recovered  the  Linear  Quadratic  Regulator  (LQR)  performance;  however,  both  seemed  to  have  a 
performance  limit.  During  the  effort  to  find  a  more  direct  measure  of  drag  reduction  with  periodic 
optimization,  we  found  that  upstream  traveling  waves  could  drive  skin-friction  drag  below  laminar. 
LQG  controllers  for  Rayleigh-Benard  convection  suppression  were  more  successful;  reduced-order 
controllers  successfully  stabilized  convection  in  linear  and  non-linear  simulations  of  a  heated  layer 
of  fluid. 


2.1  LQG/Loop  Transfer  Recovery  Matches  LQR  in  Flow  Control 
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Figure  1:  Drag  evolution  of  LQG  and  LQR  controlled  turbulent  flow.  7  is  the 
regulator  cost  function  tuning  parameter,  which  weights  the  control  input  against 
the  cost  function  of  the  perturbed  streamwise  wall-shear  stress. 

The  Navier-Stokes  equations,  linearized  about  a  laminar  Poiseuille  flow  profile,  are  used  as  a 
basis  for  the  design  of  reduced-order  controllers.  The  decomposition,  a  two-dimensional  Fourier 
expansion  in  the  streamwise  and  spanwise  directions  and  a  Galerkin  projection  in  the  wall-normal 
direction,  decoupled  the  dynamical  equations  by  each  Fourier  wavenumber.  LQG  synthesis  is 
applied  to  a  few  of  these  independent  sub-systems,  in  effect  applying  control  to  suppress  pertur¬ 
bations  only  at  certain  Fourier  wavenumbers.  Furthermore,  each  sub-system  is  further  reduced 
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by  selecting  only  the  most  controllable  and  observable  modes.  Thus,  the  controller  equations  are 
reduced  98%  in  size.  Application  on  fully  nonlinear  turbulent  DNS  of  a  channel  flow  have  been 
moderately  successful  with  approximately  15%-18%  reduction  in  drag  [4].  Linear  mechanisms  are 
demonstrated  to  be  an  important  part  of  sustaining  turbulence  [9]. 

Comparison  of  the  un-rotated,  full-state  linearized  equations  evolution  and  nonlinear  DNS  with 
sufficiently  small  amplitudes  revealed  slight  differences  in  the  truncation  errors  present  in  both  sim¬ 
ulations.  The  differences  were  significant  enough  to  prevent  the  estimation  filter  from  accurately 
reconstructing  the  plant  state.  The  estimation  was  significantly  improved  by  accounting  for  the 
differences  in  truncation  errors,  and  the  LQG  controllers  were  able  to  perform  about  as  well  as 
the  full  state  feedback,  LQR  controllers  (see  Figure  1).  However,  in  experiments  with  different 
cost  functions,  such  as  minimal  perturbation  energy  and  reducing  the  coupling  term  to  the  Squire 
equations,  the  resulting  controllers  were  unable  to  gain  any  further  reductions  in  drag. 


2.2  Upstream  Travelling  Wave:  sub-laminar  drag 

Our  linear  control  design  was  based  on  the  notion  that  if  the  perturbations  can  be  removed  at  all 
significant  wavenumbers,  the  the  indirect  effect  would  be  a  reduction  in  mean  viscous  drag.  Since  it 
became  clear  that  this  indirect  approach  to  drag  reduction  would  be  of  limited  effectiveness,  a  more 
direct  approach  that  links  perturbations  in  each  wavenumber  pair  to  drag  through  the  nonlineari¬ 
ties  in  the  the  governing  equations  was  investigated.  Noting  that  our  LQR  controllers  transiently 
induced  below  laminar  drag  [10]  with  certain  initial  flow  fields,  we  attempted  to  determine  if  there 
were  a  periodic  process  in  which  average  drag  could  be  sustained  below  that  of  laminar  flow.  Initial 
attempts  using  a  sequence  of  two  controllers,  one  that  induced  below  laminar  drag  and  the  other  to 
recover  the  initial  flow  conditions,  did  not  result  in  reduced  averaged  sustained  drag  below  laminar 
flow.  Therefore,  we  took  a  direct  approach  by  developing  a  numerical  code  to  determine  the  opti¬ 
mal  periodic  open-loop  control.  Using  our  DNS  code  as  the  plant,  perturbations  were  made  in  the 
parameterized  controls  as  well  as  the  initial  flow  conditions  for  a  small  set  of  wave  number  pairs. 
Numerical  gradients  were  constructed  and  became  inputs  to  an  accelerated  gradient  optimization 
technique. 

It  might  appear  that  this  effort  would  not  produce  sustained  drag  below  laminar  flow  based  on 
a  conjecture  proposed  by  Bewley  and  Aamo  [11]: 

Conjecture:  The  lowest  sustainable  drag  of  an  incompressible  constant  mass-flux 
channel  flow,  when  controlled  via  a  distribution  of  zero-net  mass-flux  blowing/suction 
over  the  no-slip  channel  walls,  is  exactly  that  of  the  laminar  flow. 

However,  in  the  process  of  simplifying  this  complex  numerical  scheme,  we  discovered  that  two- 
dimensional  waves  of  wall-bounded,  wall-normal  blowing  and  suction  made  to  travel  upstream  can 
reduce  skin-friction  drag  on  a  sustained  basis.  This  traveling  wave  control  is  defined  in  physical 
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space  as 


vw(x,  t)  —  a  cos(kx(x  —  c£))  , 

where  a  is  a  scalar,  c  is  the  convective  speed  normalized  to  the  mean  centerline  velocity,  kx  is  the 
streamwise  wavenumber,  vw  denotes  the  wall-normal  velocity  at  the  wall,  i.e.  wall-normal  suction 
or  blowing,  and  is  applied  to  both  the  upper  and  lower  walls  of  the  channel  in  varicose  mode,  where 
suction  (and  blowing)  occurs  at  the  same  streamwise  location  at  both  walls.  By  convention,  the 
control  on  the  upper  wall  will  have  a  minus  sign  to  denote  blowing  (positive  sign  for  suction)  and 
the  opposite  for  the  lower  wall. 


Figure  2:  Time  history  of  drag  in  DNS  at  Re=2000.  Upstream  traveling  wave 
control  applied  at  kx  =  0.5,  speed  —2  (upstream).  Initially  turbulent  simulations 
are  full  three-dimensional,  while  initially  laminar  simulations  are  two-dimensional 
to  reduce  compuational  requirements. 


Although  the  exact  mechanism  is  not  fully  understood,  upstream  traveling  waves  have  been 
found  to  be  able  to  force  sustained  below  laminar  drag  in  simulations  initialized  by  both  laminar 
and  turbulent  flows  (Figure  2).  The  efficiency  of  these  controllers  as  the  ratio  of  the  power  saved 
minus  the  input  power  to  input  power,  r],  is 

Powersaved  -  Power input 

V  = - 5 - —  ,  (1) 

Power  input 
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Here  Powersaved  is  the  drag  reduction  from  the  laminar  solution,  and  Powerinput  is  calculated  from 
an  energy  balance  on  the  channel  volume  and  represents  the  amount  of  power  needed  to  push  or 
draw  fluid  out  of  the  channel. 

Preliminarily  calculations  of  efficiency  for  below  laminar  drag  initialized  from  laminar  flow  is 
-q  =  —0.50  for  the  0.1  amplitude  case  in  Figure  2  and  //  =  —0.54  for  the  other.  Initialized  from  the 
turbulent  channel  flow,  the  efficiencies  improves  to  more  favorable  values,  q  =  4.67  (amplitude  = 
0.1)  and  q  =  7.46  (amplitude  =  0.25).  Further  research  is  required  to  improve  the  efficiency. 

The  recent  derivation  by  Fukagata  et  al.  [12,  13]  provides  some  explanation  about  the  mech¬ 
anism  which  sustains  drag  reduction  to  below  laminar  levels.  For  a  channel  flow,  they  derived  an 
explicit  relationship  between  skin-friction  drag  value  and  the  Reynolds  shear  stress.  The  following 
form  follows  from  their  derivation: 

Q-Q  f*  1 

{D)oo  =  (D)i  H — —  J  uvydy,  (2) 

where  (D) ^  represents  the  total  skin -friction  drag,  (D)i  denotes  the  laminar  drag  value,  and  Re 
is  the  Reynolds  number  of  flow  centerline  velocity,  u  is  the  streamwise  perturbation  velocity,  v 
is  the  wall-normal  perturbation  velocity,  uv  denotes  the  Reynolds  shear  stress,  and  y  denotes  the 
wall-normal  direction.  While  Fukagata  et  al.  have  reported  a  conceptual  feedback  controller  using 
Lorentz  body  forces  to  achieve  sub-laminar  skin -friction  drag  [13],  our  open-loop  control  drives 
the  integral  of  the  Reynolds  shear  stress  to  become  negative  with  more  realizable  control  inputs, 
i.e.  surface  blowing  and  suction. 

In  order  to  understand  how  the  control  is  able  to  drive  the  Reynolds  shear  stress  negative,  the 
linear  mechanism  associated  with  this  phenomena  is  exploited.  Solutions  to  the  Orr-Sommerfeld 
equation  of  the  channel  flow  reveal  negative  distributions  of  the  Reynolds  shear  stress  when  sub¬ 
jected  to  upstream  traveling  waves.  Using  the  linearized  equations,  we  can  easily  find  the  expected 
drag  reductions,  i.e.  J^uvydy,  for  any  wave  speed  (Figure  3)  at  steady-state.  We  can  see  that 
while  the  channel  responds  very  adversely  to  downstream  traveling  waves,  at  certain  upstream 
speeds,  i.e.  c  <  0,  the  drag  will  be  reduced.  Furthermore,  it  suggests  that  there  is  some  opti¬ 
mal  speed  at  which  the  most  amount  of  drag  reduction  can  be  achieved.  The  peak  that  appears 
at  downstream  wave  speeds  (approximately  c  =  0.515  for  wave  at  kx  =  0.5)  corresponds  to  the 
most  observable  and  controllable  mode.  While  adverse  for  drag,  the  downstream  traveling  waves 
are  also  worth  studying  for  desirable  effects  in  certain  applications.  For  example,  it  appears  that 
an  optimized  downstream  traveling  wave  can  be  used  to  prevent  or  delay  separation  in  turbulent 
boundary  layers  subject  to  a  strong  adverse  pressure  gradient  (e.g.  diffusers  or  flow  over  airfoil). 

As  shown  in  Figure  4,  for  small  amplitudes  the  linear  equations  predict  very  well  the  amount 
of  drag  reduction  that  will  be  achieved  in  the  nonlinear  DNS.  Thus  we  believe  that  the  linearized 
equations  provide  a  good  approximation  in  which  to  analyze  the  channel  flow  response  to  travel¬ 
ing  waves.  This  provides  a  computationally  inexpensive  avenue  for  further  analysis  and  control 
development  for  improved  drag  reduction.  Details  of  our  preliminary  analysis  are  given  in  [14] 
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Figure  3:  The  weighted  Reynolds  shear  stress  integral  as  a  function  of  the  control 
wave  speed  (normalized  by  the  centerline  velocity)  for  a  Poiseuille  channel  flow  at 
Re  =  1500. 


Figure  4:  Comparison  of  drag  reduction  predicted  by  linearized  equations  and  found 
in  nonlinear  DNS. 

(Appendix  D). 


2.3  Robust  Feedback  Control  of  Rayleigh-Benard  Convection 

Over  the  last  three  years  we  have  developed  feedback  control  methods  for  the  suppression  of 
thermal  convection,  that  have  potential  applications  to  material  processing,  molding  and  crystal 
growth.  These  controllers  have  been  shown  to  be  quite  robust  within  very  realistic  simulations. 
These  controllers  have  matured  sufficiently  to  warrant  experimental  verification. 
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Although  a  proportional  feedback  controller  has  shown  certain  effectiveness  in  increasing  the 
Rayleigh  number  at  which  convection  starts,  much  better  performance  is  obtainable  by  using  a 
LQG  controller.  The  LQG  provides  a  systematic  and  optimal  way  to  synthesize  a  high-order  com¬ 
pensator,  that  guarantees  a  high  stability  margin.  Comparable  performance  by  a  PID  method  will 
require  strategic  placements  of  compensator  poles  and  zeroes  for  a  multiple-input  and  multiple- 
output  system,  which  is  a  difficult  task  if  the  system  has  large  order. 

A  closed-loop  linear  stability  analysis  was  published  in  [5].  Detailed  stability  boundary  curves 
are  computed.  Since  the  study  was  to  show  the  controller  performance,  only  the  Prandtl  number, 
Pr  =  7.0,  was  considered.  This  value  corresponds  to  water  at  room  temperature  as  the  working 
fluid.  The  sensor-actuator  model  consists  of  three  embedded  sensor  planes  and  one  actuator  plane 
coincided  with  the  lower  wall.  Both  sensing  and  actuating  are  assumed  to  be  spatially  and  tem¬ 
porally  continuous.  In  this  configuration  the  closed-loop  stability  curves  indicate  that  convection 
suppression  up  to  Rayleigh  number  (. Ra )  equal  to  14.5  times  the  critical  value  ( Rac0  =  1707.76)  is 
achievable.  However,  these  are  linear  stability  results  and  require  infinitesimally  small  amplitudes 
of  convection  in  the  fluid  layer.  In  applications,  this  requirement  is  probably  not  achievable. 

To  relax  the  restriction,  we  replace  the  linear  plant  model  with  a  three-dimensional  nonlin¬ 
ear  model  [6]  (Appendix  A).  Developed  using  the  pseudo-spectral  method,  this  model  provides 
a  realistic  simulation  of  a  heated  fluid  layer.  The  nonlinear  time-domain  simulations  show  that 
a  finite-amplitude,  steady-state  convection  in  a  layer  of  fluid  at  Ra  =  104  can  be  removed  by 
applying  the  same  linear  optimal  controller  designed  in  the  previous  study  [5]. 

In  [7]  (Appendix  B),  two  important  advances  are  made:  (i)  a  reduced-order  compensator  is  de¬ 
veloped  and  (ii)  a  gain-schedule  approach  is  proven  in  concept  by  simulation.  For  order  reduction, 
results  show  that  at  least  8  Chebyshev  modes  are  required  in  the  compensator  for  effective  perfor¬ 
mance.  By  using  the  gain  scheduling  approach  which  involves  a  two-step  Ra  profile,  we  are  able 
to  raise  the  suppressed  no-motion  state  to  about  12  times  Rac0.  The  closed-loop  linear  stability 
limit  occurs  at  about  14.5  times  R.a,,o,  which  puts  an  upper  bound  to  the  nonlinear  simulation.  In 
[8]  (Appendix  C),  the  models  have  been  improved  by  eliminating  the  idealized  thermal  boundary 
conditions  and  replacing  them  by  walls  of  finite  thickness.  As  for  deviational  errors,  the  compen¬ 
sator  model  is  run  at  the  nominal  values,  while  the  corresponding  parameters  in  the  plant  model  are 
allowed  to  change  from  their  nominal  values  from  which  stability  margins  are  computed.  Stability 
margins  are  also  determined  for  the  Rayleigh  number,  the  wavenumber,  the  actuator  delay,  and  the 
sensor-plane  depth  uncertainty.  The  conservatism  of  these  robust  stability  margins  are  investigated 
via  a  fully  nonlinear,  3D  simulation.  As  the  conclusion  of  the  study,  we  believe  the  concept  can 
now  be  implemented  as  laboratory  experiments. 
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3  Conclusions  and  Future  Work 


It  is  clear  from  the  results  of  this  study  that  Linear  Quadratic  control  laws  can  be  effective  in 
controling  nonlinear  phenomenon.  However,  there  is  a  limit  to  its  performanace.  While  the  control 
can  be  highly  effective  if  a  direct  measure  of  the  desired  quantities  is  available  such  as  with  the 
Rayleigh-Benard  convection  suppression  problem,  the  linear  controllers  are  severely  limited  when 
the  desired  objective  is  only  available  indirectly,  such  as  drag  reduction  in  channel  flows. 

The  effect  of  traveling  waves  of  wall-bounded,  wall-normal  suction  and  blowing  definitely 
warrants  further  study.  While  using  the  same  actuation  that  appeared  limited  with  the  LQR/LQG 
controllers,  a  simple  open-loop  scheme  is  able  to  demonstrate  dramatic  drag  reduction,  even  driv¬ 
ing  the  drag  below  sub-laminar  levels.  The  results  found  in  this  study  are  very  preliminary.  A 
better  understanding  of  the  phenomonen  as  well  as  the  flow  field  induced  by  these  traveling  waves 
will  be  required  before  improvements  can  be  made  to  the  control  schemes,  both  in  terms  of  further 
reduction,  flow  stability,  and  efficiency  of  control. 
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We  study  by  a  fully  nonlinear  three-dimensional  pseudospectral  time-splitting 
simulation,  the  feedback  control  of  a  layer  of  fluid  heated  from  below.  The  initial 
condition  corresponds  to  a  steady  large-amplitude  preferred  convection  state  obtained 
at  Prandtl  number  of  7.0  and  Rayleigh  number  of  104,  which  is  about  six  times  the 
Rayleigh  critical  value.  A  robust  controller  based  on  the  LQG  (linear-quadratic- 
Gaussian)  synthesis  method  is  used.  Both  sensors  and  actuator  are  thermal-based, 
planar,  and  assumed  to  be  continuously  distributed.  The  simulated  results  show  that 
large-amplitude  steady-state  convection  rolls  can  be  suppressed  by  the  linear  LQG 
controller  action.  The  Green’s  function  of  the  controller  gives  the  shape  of  the  control 
action  corresponding  to  a  point  measurement.  In  addition,  for  Rayleigh  numbers 
below  the  proportional  feedback  control  stability  limit,  this  controller  appeared  to 
be  effective  in  damping  out  steady-state  convection  rolls  as  well.  However,  in  a 
region  very  near  the  proportional  control  stability  limit,  proportional  control  action 
induces  subcritical  g-type  hexagonal  convection,  which  is  obtained  here  through  direct 
simulations.  Note  that  well  above  this  proportional  control  limit,  the  LQG  still  damps 
out  all  convection.  The  nonlinear  plant  model  is  validated  by  comparing  check  cases 
with  published  results. 


1.  Introduction 

Active  suppression  of  onset  of  convection  in  a  layer  of  fluid  has  potentially 
important  applications  in  improving  the  material  that  goes  through  solidification  in 
a  mould.  For  instance,  during  the  growth  phase  of  large  silicon  wafers  or  composite 
materials,  a  large  thermal  gradient  typically  causes  undesirable  convective  motions 
in  the  melt.  To  understand  the  active  control  of  the  realistic  manufacturing  process, 
an  idealized  system  is  an  important  starting  point.  To  this  end  Rayleigh-Benard 
convection  (RBC)  is  ideal  for  vigorous  theoretical  analysis. 

Many  theoretical  studies  have  employed  the  linear  feedback  control  to  increase 
the  stability  threshold  of  the  purely  heat  conductive  state  so  that  no  convection 
occurs  despite  the  presence  of  a  large  thermal  gradient  (Tang  &  Bau  1993,  1994, 
1998a,  b\  Howie  1997a,  h,  c,  2000;  Or,  Cortelezzi  &  Speyer  2001).  These  studies  used 
the  linear  plant  model  and  employed  a  simple  controller  using  the  proportional 
feedback.  The  implantable  sensor  and  actuator  are  assumed  to  be  of  the  thermal 
type  and  continuously  distributed  spatially  on  the  horizontal  plane.  Analysis  as 
well  as  experimental  results  in  general  indicate  that  the  proportional  controller  will 
stabilize  the  basic  state  up  to  Rayleigh  number  (Ra)  of  3  to  4  times  its  critical  value 
of  the  basic  state  (see  Tang  &  Bau  1994;  Howie  1997a).  Furthermore,  as  shown  in 
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Tang  &  Bau  (1994),  a  controller-induced  oscillatory  instability  occurs  at  a  large  gain.  A 
linear-quadratic-Gaussian  (LQG)  controller  has  also  been  studied  (Or  et  al.  2001)  to 
increase  the  region  of  stabilization  and  with  a  higher  margin  of  robustness.  First,  the 
stability  limit  can  be  raised  to  about  14  times  the  critical  value  of  Ra.  Secondly, 
the  gain  and  phase  margins  about  the  design  point  of  the  controller  appear  adequate 
for  practical  implementation. 

To  develop  a  control  design  to  be  implementable  for  applicational  processes  (such 
as  for  crystal  growth  or  a  melt),  it  is  crucial  to  understand  the  control  process  for 
simpler  geometry  and  material  properties.  We  have  been  focused  on  an  Oberbeck- 
Boussinesq  model  for  a  horizontal  layer  of  fluid.  The  plant  dynamics  is  known  as 
Rayleigh-Benard  convection  (RBC)  (Cross  &  Hohenberg  1993).  As  a  first  step,  the 
performance  of  the  linear  controller  design  for  the  linear  plant  dynamics  is  reported  in 
Or  et  al.  (2001).  In  this  paper,  as  a  step  further,  the  focus  is  turned  to  the  performance 
of  the  linear  controller  design  for  the  fully  nonlinear  plant  dynamics. 

It  is  well  known  that  in  a  large  layer  of  heated  fluid,  convection  occurs  as  a 
steady  pattern  of  two-dimensional  rolls.  The  two-dimensional  convection  rolls  and 
the  stability  properties  were  investigated  in  detail  by  Clever  &  Busse  (1974)  and 
Busse  &  Clever  (1979).  For  the  heated  layer  corresponding  to  Ra>Rac0  ( Rac0  is 
computed  theoretically  to  have  the  value  of  1707.762  up  to  3  decimal  places),  the 
stable  roll  pattern  occurs  only  within  a  band  of  wavenumber  centred  approximately 
about  the  critical  wavenumber  ac  =  3.117.  Within  the  stable  band  the  rolls  realized 
do  not  necessarily  have  a  preferred  length  scale.  Indeed,  their  wavelength  appears 
to  be  dictated  by  the  initial  conditions  used  to  select  the  rolls  and  by  the  manner 
that  the  basic  state  temperature  is  prescribed  spatially  and  temporally.  The  band  is 
bounded  on  both  sides  by  instabilities  that  pertain  to  changing  the  wavelength  of 
the  rolls  but  not  changing  the  planform.  As  the  induced  rolls  acquire  a  wavelength 
too  large  or  too  small,  an  instability  will  occur  to  shift  their  length  scale  back 
to  a  value  close  to  the  critical  value.  As  the  value  of  Ra  increases,  the  rolls  will 
at  some  point  becomes  unstable  and  the  convection  structure  will  converge  to  a 
pattern  with  more  complex  spatial  or  temporal  structure.  The  exact  value  of  Ra  at 
which  the  transition  occurs  is  wavenumber  dependent.  For  a  Prandtl  number  (Pr)  of 
7.0,  for  instance,  the  two-dimensional  rolls  become  unstable  to  a  three-dimensional 
bimodal  convection  at  roughly  Ra  «  25 Rac0  at  the  wavenumber  of  about  2.0  (see  the 
experimental  observations  presented  in  figure  11,  Busse  &  Clever  1979).  The  transition 
highlights  a  sufficiently  strong  thermal  boundary-layer  effect,  made  possible  at  Ra 
values.  The  transition  to  three-dimensional  convection  occurs  at  a  significantly  higher 
Ra  threshold  than  the  closed-loop  stability  limit  of  14.5Pur0  based  on  the  linear  LQG 
controller  (Or  et  al.  2001).  For  our  control  analysis  here,  therefore,  we  need  only  to 
consider  the  two-dimensional  rolls  as  the  initial  state  of  convection  to  be  controlled. 

Our  present  control  problem  can  be  investigated  most  effectively  by  the  use  of  time- 
domain  analysis.  A  three-dimensional  fully  nonlinear  pseudospectral  model  using  a 
time-splitting  integration  method  is  developed,  based  on  the  Oberbeck-Boussinesq 
equations.  The  proportional  feedback  controller  is  easily  implementable.  This  case 
provides  the  check  cases  for  code  validation  purposes.  Certain  flow  patterns  that 
are  known  to  be  induced  by  the  controller  effects,  such  as  the  oscillation  mode 
(Tang  &  Bau  1994)  and  the  g-type  hexagons  (Or  &  Kelly  2001),  can  be  obtained 
here  from  the  direct  numerical  simulations  and  compared  with  those  reported  from 
previous  analyses.  In  §2,  the  nonlinear  plant  model  and  the  LQG  controller  will 
be  briefly  described.  The  results  will  be  presented  in  §3,  followed  by  the  conclusion 
in  §4. 
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2.  Mathematical  formulation 

2.1.  Nonlinear  plant  model  and  numerical  solution 
The  nonlinear  plant  model  is  governed  by  the  Oberbeck-Boussinesq  equations  for  a 
horizontal  layer  of  fluid.  In  the  non-dimensional  form  d,  d2 /at,  </d,  t</d2,  p(K/d )2  and 
AT  are  used  as  the  scales  of  length,  time,  velocity,  vorticity,  pressure  and  temperature, 
where  d  is  the  layer  thickness,  k  and  p  are  the  mean  thermal  diffusivity  and  density 
of  the  fluid,  and  AT  is  the  temperature  difference  between  the  upper  and  lower  wall 
in  the  purely  conductive  basic  state.  The  governing  non-dimensional  equations  are, 

PN1  d,v  =  Pr~lv  x  (o  +  kRaO  —  Vji  +  V2r,  (2.1) 

3,0  =  -v  ■  V6>  +  w  +  V26>,  (2.2) 

V  •  v  =  0,  (2.3) 

where  v  =  (u,  v,  w)  is  the  velocity  vector  field,  co  =  Vx  v  is  the  vorticity,  n  =  p+vv/2 
is  the  pressure  head,  0  is  the  perturbation  temperature  and  k  is  the  unit  vector  in  the 
z-di  recti  on.  The  two  external  parameters  are  Rayleigh  and  Prandtl  numbers,  given 
by  Ra  =  gATd2/vK  and  Pr  =  v/k  where  v  is  the  mean  kinematic  viscosity.  The 
continuity  equation  (2.3)  applies  only  when  the  flow  is  incompressible. 

The  velocity  field  is  assumed  to  be  non-permeable  and  non-slip  at  the  upper  and 
lower  walls,  thus  subject  to 

v(x,  y,  0,  t)  =  0,  v(x,  y,  1,  t)  =  0.  (2.4) 

The  temperature  field,  on  the  other  hand,  is  assumed  to  satisfy  the  isothermal  condition 
at  the  upper  wall.  The  lower  wall  is  non-isothermal  owing  to  the  action  of  the  thermal 
actuation.  It  is  assumed  that  a  control  temperature  9c(x,y,t )  can  be  imposed.  The 
upper  and  lower  thermal  boundary  conditions  for  the  perturbation  field  are  therefore 

0(x,  y,  1,  t)  =  0,  9(x,  y,  0,  t)  =  9c{x,  y,  t).  (2.5) 

In  order  to  perform  the  feedback  control,  the  perturbation  temperature  field  has  to 
be  measured  in  the  fluid.  In  our  model,  three  sensor  planes  are  embedded  in  the 
layer  at  carefully  chosen  levels  at  z  =  zs  (with  s  =  1,2  and  3).  For  analysis  purposes, 
these  sensor  planes  are  assumed  to  exert  no  blockage  effects  on  the  flow  field.  They 
measure  the  planar  temperature  distribution  in  the  layer, 

9{x,  y,  zs,  t)  =  9s{x,  y,  t),  s  =  1,2,3.  (2.6) 

Assuming  a  continuous-distributed  sensor,  9s(x,  y,t )  are  known  at  sampled  points 
and  time. 

In  the  numerical  scheme,  the  dependent  variables  u,  v,  w,  p  and  9  are  expressed  by 
the  following  truncated,  triple  sums, 

^  (  Mkmn 

V  I  N  ^  K  ^  Vkmn 

w  (x,  y,  z,  t)  =  Re  EE  E  U’kmn  ( t)Tn(z )  exp (i{kaxx  +mayy ))  >  (2.7) 

P  n= 0  k= 0  m=— M+l  Pkmn 

_  @  _  <  _  @kmn  _  > 

where  Re  denotes  the  real  part  of  the  sum,  ax  and  ay  are  the  fundamental 
wavenumbers  in  the  x-  and  y-directions,  respectively.  The  asymmetric  treatment 
of  the  indices  k  and  m  reduces  the  number  of  coefficients  by  half  because  the  velocity, 
pressure  and  temperature  are  real  dependent  variables  (see  Marcus  1984).  These  two 
parameters  are  prescribed  in  the  model.  The  functions  T„(z)  («  =  0,1,...)  denote  the 
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Chebyshev  polynomials.  Note  that  a  linear  coordinate  transformation  is  implicitly 
assumed  to  convert  the  Chebyshev  function  domain  between  +1  and  —1  to  our 
physical  range  0  ^  z  <  1  •  The  actuator  and  sensor  temperatures,  6C  and  9S  (s  =  1,2 
and  3),  are  planar  (two-dimensional).  They  are  expanded  in  double  series  in  a  similar 
fashion, 

%(z i,  t) 

9S(Z2,  t) 

0S(Z3,  t) 

_0C(0 ,0. 

In  our  terminology,  the  lower  thermal  boundary  condition  (2.5)  and  the  sensor 
equations  (2.6)  are,  respectively,  the  input  to  and  output  from  the  nonlinear  plant 
model. 

The  nonlinear  equations,  together  with  the  boundary  and  the  continuity  equations 
are  then  solved  numerically  by  using  the  pseudospectral,  time-splitting  integration 
technique  (Gottlieb  &  Orszag  1977;  Canuto  et  al.  1986;  Bodenschatz,  Pesch  &  Ahlers 
2000).  Marcus  (1984)  provided  a  detailed  description  of  the  implementation  of  the 
method  for  the  Taylor-vortex  flow  simulations.  Using  the  time-splitting  method,  an 
integration  time  step  is  split  into  three  fractional  steps.  The  first  is  a  nonlinear 
fractional  step,  typically  done  using  an  explicit,  second-order  Adams-Bashforth 
scheme, 

j,iV+l/3  =  VN  +  a,3| yv  x  +  PrRa9Nk ]  _  At  1  iyv-1  x  MN- 1  +  prRa9N-ik]t  (2.9) 

eN+ 1/3  =eN  -  Atl[v" •  -  wN )  +  a t\[vN~x  •  ve"-1  -  u,"-1].  (2.10) 

The  superscript  N  here  denotes  the  time  step  and  is  not  to  be  confused  with 
the  truncation  number  for  the  vertical  dependence.  A  significant  fraction  of  the 
total  computation  load  occurs  in  computing  the  nonlinear  terms.  In  the  collocation 
space,  the  nonlinear  terms  are  computed  spatially  by  point-by-point  multiplications. 
However,  fast  Fourier  transforms  (FFT)  and  inverse  fast  Fourier  transforms  (IFFT) 
have  to  be  used  to  convert  the  field  back  and  forth  between  the  collocation  and  the 
Chebyshev-Fourier  spaces.  The  fast  Fourier  transform  (FFT)  and  inverse  fast  Fourier 
transform  (IFFT)  routines  are  obtained  from  the  library  of  the  Numerical  Recipes 
(Press  et  al.  1992),  with  some  minor  modifications.  For  validation,  these  routines  have 
been  checked  against  the  standard  Matlab  FFT  and  IFFT  functions  and  match  up 
to  15  decimal  places.  For  typical  flow  fields,  the  truncation  errors  from  FFT  and 
IFFT  due  to  aliasing  are  mainly  small  (Marcus  1984;  Press  et  al.  1992).  We  note 
that,  however,  that  the  FFT  method  can  still  be  computationally  demanding  for  high- 
resolution  solutions.  The  pseudospectral  method  is  generally  known  to  be  efficient. 
There  also  exist  other  efficient  methods  not  using  the  transforms,  for  instance,  the 
reduced-order  Galerkin  method  (Howie  1996). 

After  obtaining  the  ( N  +  l/3)th  fractional  step  with  the  Adams-Basforth  scheme, 
we  compute  the  ( N  +  l)th  step  from  the  following  equation, 

(1  -  PrAtV2)vN+1  =  uw+1/3  -  PrAtVn,  (2.11) 

subject  to  V  •  vJV+1  =  0.  It  is  noted  that,  in  general,  V  •  uJV+1/3  =f=  0.  The  most 
straightforward  procedure  for  solving  (2.11)  appears  to  be  splitting  the  equation 
into  a  pressure  step  and  a  viscous  step  (we  refer  to  it  as  the  direct  approach).  In 
the  pressure  step,  the  flow  field  subject  to  a  no  normal-flow  boundary  condition  at 
the  walls  can  be  solved  from  a  Poisson  equation,  based  on  the  property  that  the 
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pressure  field  is  irrotational  and  the  flow  field  satisfies  the  continuity  constraint  (2.3). 
Next,  a  diffusive  fractional  step  completes  the  solution  of  the  fractional  velocity  and 
temperature  fields  by  prescribing  the  no-slip  and  thermal  boundary  conditions  at  the 
walls.  As  simple  as  it  appeared,  the  scheme  had  problems  computing  the  correct  flow 
field.  In  his  numerical  simulation  of  Taylor  vortex  flow,  Marcus  (1984)  reported  large 
boundary  errors  using  this  direct  approach.  In  his  discussions  it  was  argued  that  the 
shear  may  play  a  role  and  it  is  not  clear  whether  a  similar  problem  will  occur  for 
RBC.  In  our  study  we  have  applied  the  direct  scheme  in  our  preliminary  simulations 
and  observed  large  errors  even  for  the  open-loop  simulations.  Thus,  it  appears  that 
the  problem  is  common  to  both  Taylor  vortex  flows  and  RBC.  For  more  detail  about 
the  cause  of  the  large  boundary  errors  in  the  direct  approach,  we  refer  to  Marcus 
(1984).  Marcus  identified  the  source  of  errors  and  developed  a  procedure  to  correct 
it.  His  remedy  is  to  further  split  the  fractional  solution  into  a  complementary  and  a 
particular  solution  so  that  the  boundary  conditions  and  the  continuity  equations  are 
satisfied  numerically.  The  procedure,  however,  involves  the  additional  computation 
of  several  Green’s  functions  and  seems  elaborate.  Since  the  boundary-value  problem 
corresponding  to  (2.11)  is  linear,  we  anticipate  there  are  simpler  alternative  approaches 
to  resolve  the  numerical  difficulty.  Here,  we  solve  the  problem  involving  the  pressure 
and  viscous  forces  as  a  single  step,  without  splitting  the  pressure  and  viscous  terms. 
First,  we  use  the  continuity  equation  as  the  constraint  and  eliminate  the  two  horizontal 
velocity  components  in  favour  of  the  vertical  component.  Secondly,  we  obtain  the 
solution  of  the  boundary-value  problem  for  w  and  6.  Finally,  we  recover  u  and  v  from 
w,  Fourier  mode  by  Fourier  mode,  again  using  the  continuity  equation.  This  scheme 
seems  significantly  simpler  and  has  been  tested  here  to  be  effective.  Because  of  the 
simplicity,  it  is  worth  the  description  as  an  alternative  approach  to  the  time-splitting 
procedure. 

By  eliminating  pressure  from  (2.11),  we  obtain  a  single  scalar  equation  governing  w, 

(1  -  PrAtV2)V2dzwN+1  =  Vi dzwN+1/3  -  d2z(dxuN+l/3  +  dyvN+1/3).  (2.12) 

The  equation  above  is  integrated  in  z,  this  gives 

(1  -  PrAtV2)V2wN+1  =  V2±wN+l/ 3  -  (d2zuN+1/3  +  d2zvN+1/ 3).  (2.13) 

The  integration  constant  is  zero  because  of  the  non-slip  boundary  condition.  (This 
constant  will  depend  on  the  initial  conditions  when  the  case  of  free-slip  boundary 
conditions  is  considered).  Equation  (2.13)  is  of  fourth-order  spatially.  It  has  to  satisfy 
four  boundary  conditions,  as  follows, 

wN+l=0,  dzwN+1=0  at  z  =  0, 1.  (2.14) 

The  fourth-order  boundary-value  problem,  (2.12),  determines  wN+l.  After  we  have 
obtained  wN+i,  the  horizontal  velocity  components  corresponding  to  wN+1  can  be 
obtained  by  inverting  the  continuity  and  Helmholtz  equations,  Fourier  mode  by 
Fourier  mode.  In  the  expansion,  wN+1  is  given  by 
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Each  pair  of  coefficients  (uk+.  w^1)  now  satisfies  a  Helmholtz  equation 
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where  akm  =  {(kax)2  +  {may)2}1/2.  The  Helmholtz  equation  together  with  the  con¬ 
tinuity  equation  allows  us  to  solve  for  in  terms  of  w knj'1,  provided  that  akm  0. 
The  condition  akm  =f=  0  can  occur  in  the  case  of  a  free-slip  wall,  but  not  in  the  case 
of  a  no-slip  wall.  We  refer  to  the  discussion  of  Cross  &  Hohenberg  (1993,  p.  970). 
The  perturbation  temperature  field,  on  the  other  hand,  is  not  constrained  to  have 
zero  mean  field.  Using  the  continuity  equation,  we  obtain  the  horizontal  velocity 
components, 
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2.2.  The  proportional  feedback  controller 

In  the  proportional  feedback  control,  a  proportional  relationship  is  constructed 
between  the  input  and  output  of  the  plant.  As  in  the  cases  studied  by  Tang  & 
Bau  (1994)  and  Or  et  al.  (2001),  only  one  sensor  plane  is  used  and  the  control  law  in 
this  case  is 

9c(0.t)  =  -Kp9s(zs,t),  (2.19) 

where  Kp  is  a  constant  gain  and  zs  is  the  vertical  height  of  the  sensor  plane.  The 
controller  is  very  simple  for  this  case. 

2.3.  The  LQG  controller 

The  theory  and  design  of  the  LQG  controller  was  described  in  Or  et  al.  (2001).  In 
brief,  the  linear  stability  equations  of  the  Fourier-decomposed  system  of  convection 
and  the  measurement  equation  are  given  in  matrix  form,  respectively,  by 

x  =  Ax  +  Bu ,  z  =  Cx ,  (2.20) 

where  the  entries  of  the  state  vector  x  are  the  Chebyshev  coefficients  of  velocity 
and  temperature  perturbations;  u  (measured  at  plane  z  =  0)  and  z  (measured  at 
planes  zi,zi  and  zi)  are,  respectively,  the  Fourier  coefficients  of  the  planar  control 
and  measured  temperatures.  Note  that  the  Fourier-decomposed  equations  correspond 
to  wavenumber  akm  and  Rayleigh  number  Ra.  The  following  modifications  to  the 
original  formulation  of  the  controller  model  (Or  et  al.  2001)  have  been  made  here: 
(i)  the  vertical  dependence  is  expanded  in  terms  of  the  Chebyshev  polynomials 
instead  of  the  beam  functions  as  the  basis  functions.  The  expansion  procedure, 
originally  based  on  the  Galerkin  method,  has  been  converted  to  the  tau  method.  In 
the  improved  numerical  procedure,  we  obtain  the  exact  condition  D  =  0,  in  contrast 
to  the  previous  condition  that  D  — >  0  only  as  N  — >  oo.  (ii)  We  no  longer  consider  the 
wavenumber  as  a  prescribed  parameter  here.  Instead,  an  individual  modal  controller 
is  developed  for  each  set  of  wavenumbers  (, kax,may ).  There  is  a  set  of  state-space 
equations  for  each  wave  vector.  In  total,  there  are  2 (K  +  1  )M  sets  of  A,  B  and  C 
matrices  to  be  processed. 

The  LQG  controller  is  comprised  of  a  Kalman  filter  and  an  optimal  regulator. 
The  Kalman  filter  equation  and  the  optimal  regulator  equation  corresponding  to  the 
state-space  equations,  (2.20),  are,  respectively, 

x  =  A*  x  +  B"  u  +  K  f(z  —  z),  z  =  C*x,  u=—Kcx , 


(2.21) 
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Figure  1.  The  LQG  control  loop  diagram. 


where  x  is  the  estimate  state  vector.  We  distinguish  the  matrices  with  asterisk 
superscripts  to  highlight  that  the  system  is  computed  at  a  nominal  (designed) 
wavenumber  and  Rayleigh  number,  (a*km,  Ra* ).  The  Kalman  gain  vector  K f  and 
the  optimal  gain  vector  Kc  are  determined  from  separate  steady-state  algebraic 
Riccati  equations.  The  Kalman  filter  is  used  here  as  a  state  observer  rather  than  as  an 
estimator  since  no  noises  are  injected  into  the  system  simulation.  The  cost  functional, 
weighting  and  filter  parameters  chosen  for  controller  design  are  described  in  detail  in 
Or  et  al.  (2001)  which  will  not  be  repeated  here.  For  robustness,  in  the  design,  the 
Kalman  filter  input  matrix  G  has  been  set  equal  to  the  control  input  matrix  B,  a  step 
known  as  the  loop  transfer  recovery  to  recover  the  full-state  feedback  performance 
of  the  optimal  regulator.  The  weights  for  the  objective  functions,  as  well  as  the  filter 
parameters  and  the  loop  transfer  recovery  are  described  in  Or  et  al.  (2001). 

It  is  worth  noting  that  the  LQG  controller  is  a  variant  of  the  controller  when 
the  disturbance  attenuation  bound  is  infinite  (see  Rhee  &  Speyer  1991).  In  Or  et  al. 
(2001),  robustness  is  demonstrated  classically  by  having  large  gain  and  phase  margins 
in  the  closed-loop  response.  Furthermore,  if  a  full  loop  transfer  recovery  is  achieved, 
the  full-state  feedback  LQ  regulator  performance  will  have  a  robustness  of  +60°  phase 
margin  and  6db  to  infinite  gain  margin.  Since  our  system  is  non-minimal  phased, 
only  partial  recovery  is  expected.  Since  large  gain  and  phase  margins  were  obtained 
for  the  linear  system,  the  performance  of  the  LQG  controller  in  terms  of  robustness 
should  not  be  expected  to  be  significantly  different  from  that  of  the  Hm  controller. 

In  figure  1,  we  show  the  three-dimensional  nonlinear  plant  model.  The  control 
input  and  measurement  output  of  the  model  are  Fourier-Chebyshev  coefficients 
rather  than  their  collocation  point  values.  On  the  other  hand,  in  the  physical  plant 
(such  as  in  laboratory  experiments)  the  input  and  output  are  physical  temperature 
distributions.  Since  the  LQG  controller  is  formulated  in  the  modal  space,  when  the 
upper  block  represents  the  plant  instead  of  the  model,  an  FFT  and  an  IFFT  have  to 
be  performed,  respectively,  at  the  input  and  output  of  the  controller.  In  our  case,  the 
LQG  controller  takes  the  measurements  from  the  three-dimensional  nonlinear  plant 
model  (Fourier  coefficients  at  sensor  planes)  as  input  and  determines  a  control  action 
(Fourier  coefficients  at  actuator  plane)  as  output.  The  estimate  state  vector  represents 
the  vertical  structure  and  the  state  matrices  A*,  B‘.  C  and  D  are  computed  in 
terms  of  the  designed  values  of  wavenumber  and  Rayleigh  number,  a *km  and  Ra*  (see 
equation  (3.7)  of  Or  et  al.  2001). 
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The  truncation  numbers  (K  =  32,  M  =  32,  N  =  32  +  1)  considered  here  are  of 
moderate  size.  It  is  still  convenient  to  compute  and  pre-store  the  steady-state  Kalman 
gain  K f  and  regulator  gain,  Kc.  However,  it  is  not  feasible  to  pre-store  the  state 
matrices  A *  for  all  the  wavenumbers.  Instead,  we  compute  A *  for  each  set  of 
wavenumbers  at  each  time  step  in  the  time  loop.  At  each  time  step,  the  three  sensor 
plane  temperatures,  in  modal  coefficients,  9kms(Zi,  t)  (i  =  1,2  and  3)  (see  (2.8)),  are 
exported  from  the  nonlinear  plant  model.  There  are  6 (K  +  1  )M  of  such  coefficients, 
corresponding  to  wavenumbers  0  to  Kax  in  the  x  dependence  and  (-M  +  l)av  to  May 
in  the  y  dependence.  These  values  are  then  fed  into  the  controller  which  consists  of 
the  Kalman  filter  and  the  regulator.  The  controller  processes  the  information  based  on 
the  measured  data  and  determines  the  control  output  in  terms  of  a  set  of  2(K  +  \)M 
modal  coefficients  for  9km,c( 0-  t).  These  values  are  then  input  into  the  nonlinear  model 
through  the  lower-wall  boundary  condition. 

2.4.  Green’s  function  for  point  sensors  and  actuators 
In  some  experimental  implementations  (Tang  &  Bau  1998n),  the  sensors  and  actuators 
are  discrete  rather  than  continuous.  For  the  low-resolution  point  sensors  and  actuators 
(typically  with  spacing  between  array  points  of  0(d)),  it  is  desirable  to  stack  the  arrays 
of  sensor  and  actuator  points  vertically  on  top  of  each  other.  Indeed,  our  result  will 
show  that  the  maximal  effect  of  actuation  caused  by  an  impulse  on  the  sensor  plane 
occurs  as  a  point  collocated  horizontally  with  the  impulse.  For  a  linear  system,  the 
controller  input-to-output  relationship  can  be  expressed  in  the  following  integral  form, 

9c(x,  y,t)  =  J  J  J  G(x,  y,  t\x',  y' ,  t')9s(x' ,  y' ,  t' )  Ax'  Ay'  At',  (2.22) 

where  9c(x,  y,  t)  and  9s(x',  y' ,  f)  are,  respectively,  the  planar  actuator  and  sensor 
temperature  fields.  Here,  (x,  y)  and  (x',  /)  denote  coordinates  for  the  actuator  and 
sensor  planes,  respectively.  The  kernel  G(x,  y,  t,  \x' ,  y',  f)  is  a  Green’s  function  (or  an 
influence  function).  The  first  three  arguments  in  G  represent  the  effect  and  the  last 
three  represent  the  cause. 

In  principle,  the  input  and  output  of  the  LQG  controller  can  be  represented 
by  a  linear  differential  operator  L.  The  precise  form  of  L  need  not  be  specified 
here,  since  for  our  purpose  the  Green’s  function  will  be  computed  spectrally.  In 
terms  of  L,  we  can  describe  some  general  properties  of  Green’s  function.  The 
input  and  output  temperatures  to  the  controller  are  governed  by  L 9C  =  9S,  subject 
to  appropriate  lateral  boundary  conditions  in  x,  y.  The  Green’s  formula  for  any 
two  arbitrary  functions  u(x,y,t)  and  v(x,y,t)  can  be  written  as  the  sum  of  an 
integral  J  f  f(uLv  —  vlA u)  Ax'  Ay'  At'  and  a  number  of  terms  evaluated  at  the  lateral 
boundaries  x  =  0,  2tt /ax  and  y  =  0,  2n/ay.  In  the  formula,  L+  is  the  adjoint  operator 
of  L.  Now  if  further  restrictions  are  imposed  on  u  and  v,  the  Green’s  formula  produces 
some  important  property  about  the  Green’s  function.  Let  u  =  G(x,  y,  t|xi,  yi,  f)  and 
v  =  G+(x,  y,  t\x2,  >’2,  t')  where  G  and  G+  satisfy,  respectively, 

LG(x,  y,  t\xi,  yu  t')  =  <5(x  —  xi)<5(y  —  yi)<5(f  —  t'),  ) 

?  (2.23) 

L+G+(x,  >\  t  |x2,  yi,  t')  =  <5(x  —  x2)<5(y  —  y2)S(t  -  t').) 

In  addition,  G  and  G+  satisfy  the  appropriate  lateral  boundary  conditions  and  adjoint 
boundary  conditions  so  that  the  boundary  terms  in  the  Green’s  formula  vanish.  The 
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Green’s  formula  becomes 


///(G+LG 


GL+G+)d.r'  dy'dt'  =  0. 


(2.24) 


Substituting  (2.23)  into  (2.24),  we  obtain  Maxwell’s  reciprocity  relationship  G 
(x2,  y2,  t\xu  yu  t')  =  G+(xi,  y2,  t\x2,  yiJ')-  In  our  problem,  the  lateral  boundary 
conditions  are  periodic.  The  differential  operators  in  x  and  y  are  even  in  dx  and 
3V.  The  linear  operator  L  is  self-adjoint,  i.e.  L  =  L+  and  the  symmetric  relationship 
holds, 

G(x 2,  y2,  t |jci,  yi,  t ')  =  G(xu  vi,  t\x2,  y2,  t’).  (2.25) 

The  symmetry  relationship  above  can  be  interpreted  as  follows :  at  a  given  time  t  >  t' , 
an  actuator  output  of  the  controller  at  (x2,  y2)  due  to  a  unit  impulsive  sensor  input 
of  the  controller  at  (x\ ,  jq)  and  time  t'  is  equal  to  the  actuator  output  at  {x2,  yi)  due 
to  a  unit  impulse  sensor  input  at  (x2,  y?)  and  time  t’ . 

Of  particular  interest  here  is  the  shape  of  the  actuator  temperature  0c(x,y,t) 
generated  by  a  unit  impulse  temperature  at  a  sensor  point  (xp,  yp),  say,  at  t  =  tp. 
The  spatial  roll-off  of  the  actuator  temperature  affects  the  spatial  resolution  of  the 
spacing  between  sensor  points.  Let  the  impulsive  measurement  be 

es(x',  /,  t')  =  fix’  -  xp)8{y’  -  yp)S(t '  -  tp ),  (2.26) 

from  (2.22)  we  obtain  the  Green’s  function 


9c(x,  y,  t)  =  G(x,  y,  t\xp,  yp,  tp). 


(2.27) 


For  each  Fourier  mode  that  corresponds  to  the  wave  vector  ( kax ,  may )  (where  —K /2  ^ 
k  ^  K /2  and  —M/2  ^  m  ^  M/2),  the  coefficient  represents  an  entry  of  measurement 
vector  z  in  the  filter  equation,  (2.21).  We  then  have 


u(z ,  kax,  moiy,  t ) 


exp  (A 


K  fC){t  —  t)z(t)  dr. 


(2.28) 


Note  that  the  homogeneous  solution  due  to  the  initial  condition  decays  rapidly  and 
does  not  contribute  for  sufficiently  large  t.  After  the  z  and  u  of  all  the  Fourier  modes 
are  computed,  a  FFT  will  transform  the  two  sets  of  coefficients  to  0s(x,  y,  zs,t)  and 
6c(x,  y,  0,  t),  respectively.  When  9s(x,y,zs,t )  is  impulsive  according  to  (2.26),  then 
(2.28)  gives  the  Green’s  function. 


3.  Numerical  results 

3.1.  Nonlinear  convection 

Above  the  value  Ra  =  Rac 0  «  1707.76,  the  no-motion  state  gives  rise  to  steady  two- 
dimensional  convection  rolls.  Depending  on  the  value  of  Pr,  these  rolls  in  turn  will 
become  unstable  at  still  higher  values  of  Ra,  making  transitions  to  two-dimensional 
oscillatory  convection  or  steady  three-dimensional  convection  depending  on  the  value 
of  the  Prandtl  number.  Cross  &  Hohenberg  (1993)  give  considerable  detail  about  the 
bifurcation  diagram. 

Before  engaging  in  the  closed-loop  numerical  simulations,  it  is  worth  performing 
some  comparisons  to  known  results,  as  check  cases  for  validating  the  nonlinear 
plant  model.  In  Clever  &  Busse  (1974),  selective  Nusselt  number  values  for  the  two- 
dimensional  convection  solution  were  published.  Table  1  shows  the  values  of  the 
Nusselt  number,  Nu,  for  several  different  values  of  Ra  at  Pr  =  0.71  and  7.0  for 
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Ra 

Pr  =  0.71 

Pr  =  7.0 

2000 

1.210  (1.212) 

1.214  (1.214) 

2500 

1.472  (1.475) 

1.475  (1.475) 

10000 

2.653  (2.661) 

2.608  (2.618) 

Table  1.  Nusselt  number  values  for  two-dimensional  rolls. 


two-dimensional  rolls  at  ax  =  3.117  (ay  =  0).  The  Nusselt  number  is  a  measure  of  the 
convective  heat  transfer,  defined  as  the  value  of  temperature  gradient  at  either  upper 
or  lower  wall. 


N 

Nu  =  1  +  0oo» 

n= 0 


d  Tn(z) 

dz 


j=0, 1 


(3.1) 


where  the  first  two  zero  indices  of  900n  correspond  to  k  =  m  =0  so  that  the  sum 
represents  the  temperature  gradient  averaged  over  the  horizontal  plane.  In  the  absence 
of  an  internal  heat  source,  the  values  of  Nu  evaluated  at  z  =  0  and  z  =  1  should  be 
equal.  Our  open-loop,  steady-state  solutions  are  obtained  at  truncation  numbers 
K  =  16,  M  =  8  and  N  =  16,  for  ax  =  3.117  and  av  =0  (transverse  rolls).  In  table  1,  the 
values  published  in  Clever  &  Busse  (1974)  are  shown  in  parentheses.  In  all  cases,  the 
difference  between  our  value  and  theirs  is  less  than  0.4%.  For  values  of  wavenumber 
ax  =  2.2  and  2.6,  respectively,  where  Pr  =  1  and  Ra  =  10000,  we  obtain  Nu  =  2.465 
and  2.548  versus  their  values  2.473  and  2.557.  We  further  note  that  Nu  should  not 
depend  on  the  orientation  of  rolls.  As  a  consistency  check,  we  compare  the  Nu  of 
our  solutions  between  the  longitudinal  (aA  =  0,  ay  =f=  0)  and  transverse  rolls  (ax  =j=  0, 
ay  =  0).  The  difference  of  the  Nu  values  is  found  to  be  less  than  0.02%. 


3.2.  Proportional  feedback  control 

We  now  turn  to  the  proportional  feedback  control  problem.  From  the  results  of 
Tang  &  Bau  (1994)  and  Or  et  al.  (2001),  oscillatory  convection  occurs  when  the 
proportional  gain  Kp  becomes  sufficiently  large.  At  Kp  =  6,  for  instance,  the  linear 
theory  at  Pr  =  7  predicts  that  an  oscillatory  instability  is  preferred  over  the  steady- 
state  rolls.  The  closed-loop  threshold  of  stability  is  ac  =  3.73  and  Ra  =  3.63 Rac0,  with 
the  frequency  of  oscillation  equal  to  20.4.  For  the  same  values  of  Ra  and  wavenumber 
we  use  the  steady  state  two-dimensional  rolls  as  the  initial  conditions  for  our  time- 
domain  simulation.  Our  results  appear  to  be  consistent  with  the  prediction  of  linear 
theory.  Figure  2  shows  the  behaviour  of  Nu  of  the  closed-loop  solutions  at  Kp  =  6 
for  two  values  of  Ra/Rac0:  at  3.55  (solid)  and  3.65  (dashed).  In  both  curves,  the 
open-loop  steady  two-dimensional  rolls  are  used  as  the  initial  condition.  These  rolls 
are  obtained  at  Ra/Rac o  =  3.65  and  ax  =  3.73  which  yield  Nu  =  2.273.  In  figure  2  the 
solid  curve  shows  stable  behaviour  whereas  the  dashed  curve  is  unstable.  The  neutral 
curve  has  Ra/Rac o  approximately  equal  to  3.60.  This  value  is  in  close  agreement 
with  the  result  of  linear  theory.  Furthermore,  the  oscillatory  behaviour  in  the  curves 
indicate  a  frequency  of  about  40.3,  again  consistent  with  eigenvalue  prediction  of 
2  x  20.4  of  the  linear  theory.  It  is  noted  that  Nu  has  a  harmonic  frequency  equal  to 
twice  the  fundamental  frequency. 

The  oscillatory  convection  appears  to  have  a  two-dimensional  roll  planform.  The 
more  interesting  finding  according  to  the  numerical  simulations  is  that  this  oscillatory 
solution  is  not  unique  for  the  given  set  of  external  parameters.  It  turns  out  when  we 
prescribe  an  additional  small  perturbation  field  in  the  y-dependence,  for  the  same 
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Figure  2.  Nusselt  number  of  oscillatory  convection. 
{a)  ( b ) 


X 


X 


Figure  3.  Closed-loop  solution:  g-type  hexagon  pattern. 

values  of  Kp,  Ra  and  ax,  the  closed-form  solution  will  not  settle  at  the  two-dimensional 
oscillatory  branch  if  the  cross-roll  perturbation  is  not  small.  For  sufficiently  large 
cross-roll  perturbations,  the  solution  will  settle  down  at  a  subcritical  branch.  In  this 
case,  the  horizontal  planform  solution  is  three-dimensional,  which  resembles  the  g- 
type  hexagons  (Or  &  Kelly  2001).  Depending  on  the  asymmetry  in  the  perturbation 
temperature,  hexagon  cells  with  sinking  motion  near  the  centre  of  the  cell  and  rising 
motion  near  the  cell  wall  is  referred  to  as  the  g-type.  For  the  f-type  hexagons  the 
opposite  is  true.  In  figure  3,  we  show  (a)  the  planform  corresponding  to  temperature 
at  the  lower  wall  (z  =  0)  and  ( b )  the  planform  corresponding  to  horizontal  velocity 
components  at  horizontal  plane  z  =  0.1  (the  velocity  components  vanish  at  the  lower 
wall  owing  to  a  non-slip  boundary  condition).  The  three-dimensional  hexagonal 
convection  is  a  steady-state  pattern  and  corresponds  to  Nu  =  1.4352.  The  hexagonal 
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Wavenumber 


Figure  4.  A  sketch  of  the  stability  boundaries  for  the  uncontrolled  layer  at  Pr  =  7.0. 


solution  induced  by  the  controller  action  has  been  studied  in  considerable  detail  (see 
Shortis  &  Hall  1996;  Or  &  Kelly  2001)  based  on  weakly  nonlinear  analysis.  Here,  we 
actually  obtain  the  solution  from  a  direct  numerical  simulation.  We  summarize  several 
important  conclusions  based  on  the  results  presented,  (i)  The  solutions  obtained  from 
our  fully  nonlinear  three-dimensional  pseudospectral  plant  model  have  been  checked 
and  agree  reasonably  well  against  known  published  results  from  other  independent 
methods,  (ii)  The  proportional  feedback  controller  induces  a  subcritical  range  of  g- 
type  hexagonal  convection,  which  appears  to  be  captured  in  the  nonlinear  simulations. 
Near  the  stability  threshold  of  the  closed-loop  system  with  sufficiently  large  gain,  both 
two-dimensional  oscillatory  convection  and  three-dimensional  steady-state  hexagonal 
convection  can  co-exist  in  the  same  parameter  region.  Next,  we  consider  the  closed- 
loop  simulation  using  the  LQG  controller. 

3.3.  Closed-loop  simulations  using  the  LQG  controller 
We  investigate  the  closed-loop  system  with  an  operating  condition  of  the  plant 
model  at  Pr  =  7.0  and  7?a  =  104.  In  the  set-up,  the  controller  gains  Kf  and  Kc  are 
steady  states,  precomputed  and  stored.  The  actual  controller  and  the  nonlinear  plant 
models  are  implemented  in  FORTRAN  and  MATLAB  .  This  controller  is  implemented 
according  to  the  description  in  §2.3. 

In  figure  4,  we  provide  a  sketch  of  the  stability  diagram  of  the  uncontrolled 
dynamics  at  Pr  =  l  (see  Busse  &  Clever  1979,  figure  1  for  the  original  plot).  The 
stability  boundary  of  the  purely  conduction  (static)  state  is  the  lowest  parabolic¬ 
shaped  curve.  At  each  Ra  above  the  minimum  of  this  neutral  curve  (supercritical), 
the  linear  theory  predicts  an  outer  band  of  wavenumbers  in  which  the  basic  state 
is  unstable.  However,  at  each  supercritical  Ra,  the  stable  finite-amplitude  convection 
occurs  in  a  narrower  band  of  wavenumbers.  At  Pr  =  7,  the  stable  finite-amplitude 
convection  in  the  inner  band  corresponds  to  steady  two-dimensional  convection  rolls. 
For  Ra  =  104  (the  dashed  line  in  figure  4),  the  inner  band  of  wavenumbers  is  bounded 
on  the  lower  side  by  the  cross-roll  instability  at  a  «  1.75  and  on  the  higher  side  by 
the  skew-varicose  instability  at  a.  «  3.5.  At  this  Ra,  the  inner  band  of  wavenumbers 
is  significantly  smaller  than  the  outer  band  obtained  from  the  linear  theory,  which 
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gives  approximately  0.74  and  9.0,  respectively.  The  stability  boundaries  are  in  general 
Pr  dependent. 

The  stable  two-dimensional  convection  rolls  are  characterized  by  a  single 
wavenumber,  but  it  can  be  any  value  within  the  inner  band.  Laboratory  experiments 
(see  Cross  &  Hohenberg  1993)  using  different  initial  conditions  had  demonstrated 
that  the  stable  pattern  can  have  a  non-unique  wavenumber.  On  the  other  hand, 
certain  experiments  performed  by  letting  Ra  vary  either  as  a  slow  function  of  time  or 
by  inducing  a  spatial  ramp  in  the  layer  thickness  indicate  that  the  rolls  are  realized 
with  a  unique  wavenumber.  Since  our  goal  here  is  to  eliminate  the  convection  pattern, 
the  detailed  properties  of  the  nonlinear  solution  do  not  concern  us  other  than  as  the 
initial  condition  for  our  closed-loop  solutions. 

The  closed-loop  simulation  is  demanding  computationally  in  the  sense  that  the 
entire  outer  band  of  wavenumbers  should  be  covered  in  the  stabilization  of  the  basic 
state.  In  our  simulation,  the  fundamental  wavenumbers  ax  and  ay  are  selected  so 
that  the  expansion  covers  the  entire  inner  band,  but  falls  short  of  the  outer  band. 
We  argue  that  this  arrangement  is  reasonable  and  we  use  the  truncation  numbers 
K  =  M  =  N  =  32.  The  nonlinearity  has  the  role  of  limiting  the  wavenumber  of  the 
convection  pattern  to  the  inner  band.  As  the  initial  condition  for  the  closed-loop 
simulation,  we  let  ax  =  1.0  in  the  open-loop  run  with  appropriate  initial  condition. 
We  obtain  a  steady,  two-dimensional  roll  pattern  with  a  wavenumber  of  3.0.  In  the 
closed-loop  simulation,  we  add  in  a  small  perturbation  of  cross-rolls  superimposed 
on  the  steady  finite-amplitude  rolls.  The  added  perturbation  assures  that  the  initial 
condition  used  is  three-dimensional. 

The  closed-loop  simulation  results  are  shown  in  figures  5(a)-5(g).  Since  the 
transition  is  two-dimensional,  it  suffices  to  reveal  the  flow  fields  by  showing  the 
cross-sectional  view  in  the  (z,  v)-plane.  In  figures  5(u)-5(c),  we  show  the  transient 
pattern  of  the  perturbation  isotherms  in  the  (x,  z)-plane  (with  basic  temperature 
subtracted).  The  three  isotherm  patterns  (figures  5a-5c)  of  the  disturbance  field  are 
snapshots  obtained  at  t  =  0,  0.05  and  0.2  diffusive  time  units,  respectively.  Note  that 
the  upper  and  lower  wall  are  both  the  perturbation  isotherm  of  zero  temperature. 
The  solid  (dashed)  lines  indicate  positive  (negative)  increments  of  temperature.  The 
same  increment  of  temperature  applies  to  all  three  panels.  Figure  5(a)  shows  the 
cross-section  of  the  steady-state  convection  rolls  used  as  the  initial  condition  at 
t  =  0.  Shortly  after  the  controller  is  turned  on  at  t  =  0,  figure  5(b)  shows  that  a 
steep  thermal  boundary-layer  pattern  develops  near  the  lower  wall  at  t  =  0.05.  This 
boundary  temperature  perturbation  possesses  an  opposite  sign  to  the  perturbation  in 
the  bulk  of  the  layer  of  fluid,  and  therefore  exerts  a  cancellation  effect,  which  tends 
to  drive  the  fluid  towards  an  isothermal  state.  Figure  5(c)  shows  at  a  later  instant 
(t  =  0.2)  that  the  isotherm  pattern  indeed  settles  towards  a  static  state.  Here,  the 
isotherms  correspond  to  a  residual  temperature  distribution  of  about  1.5%  of  the 
temperature  shown  in  figure  5(a).  The  residual  temperature  continues  to  approach 
zero  asymptotically  in  time. 

In  figures  5(d)-5(f)  we  show  the  quiver  plots  of  the  velocity  field  corresponding 
to  isotherms  in  figures  5(u)-5(c).  The  arrow  sizes  in  figures  5(e)-5(f)  are  according 
to  the  true  relative  scale.  For  illustration  of  the  flow  field  we  deliberately  magnify 
the  arrows  in  figure  5(f).  Note  that  the  velocity  rolls  are  shifted  by  a  phase  of  n/2 
relative  to  the  isotherm  rolls.  The  upward  (downward)  motion  of  fluid  is  associated 
with  the  positive  (negative)  isotherms,  as  indicated  in  the  figures.  From  figure  5(c), 
we  observe  that  the  upwelling  and  downwelling  regions  are  significantly  perturbed 
by  the  control  action.  As  a  result,  a  secondary  row  of  vortices  near  the  lower  wall 
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(a)  1.0 


z  0.5 


Figure  5.  (a)-(c)  Transient  patterns  of  isotherm  under  controller  action;  (d)-(f)  correspond¬ 
ing  velocity  quivers  of  the  flow  patterns;  (g)  Time  response  of  Nusselt  numbers  at  the  lower 
and  upper  wall  ( Ra  =  104  and  Pr  =  7.0)  (a),  ( d)t  =  0 ;(7>),  (e) 0.05; (c),  (/) 0.2. 

is  apparent.  In  figure  5(e),  the  convective  motion  becomes  so  weak  that  the  vortex 
structure  is  no  longer  visible.  Finally,  we  show  the  two  Nusselt  numbers  in  figure  5(g) 
in  time  as  the  indicator  for  convective  heat  transport.  The  lower  (solid)  and  upper 
(dashed)  curves  are  based  on  the  horizontal-mean  temperature  gradient  at  the  lower 
and  upper  walls,  respectively.  The  gradient  is  computed  normal  to  the  walls.  As  the 
thermal  actuator  action  is  switched  on,  a  large  transient  perturbation  develops  near 
the  lower  wall,  indicating  an  increase  in  local  heat  flux  from  the  actuator  action.  The 
lower  Nusselt  number  shoots  up  considerably  higher  than  the  upper  Nusselt  number 
initially  for  a  brief  duration.  Subsequently,  the  upper  Nusselt  number  is  greater  than 
the  lower  value,  as  the  heat  in  the  bulk  of  fluid  is  transferred  away.  Between  t  =  0 
and  t  =  0.474,  we  determined  through  integration  that  the  areas  under  the  curves  are 
0.5635  (solid  line)  and  0.5628  (dashed  line).  The  two  integral  values  will  converge  to 
the  same  value  in  time,  as  a  constraint  of  the  conservation  of  heat. 
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Figure  6.  (a)  Isotherms  of  the  response  of  temperature  on  the  actuator  plane  owing  to  an 
impulse  source  temperature  on  the  sensor  plane;  (b)  variation  of  the  response  temperature 
along  the  dashed  line  of  (a);  (c)  change  of  (?„„•„  with  f . 


For  a  related  drag  reduction  control  problem,  Cortelezzi  &  Speyer  (1998)  developed 
a  robust  reduced-order  controller.  It  is  beyond  the  present  scope  to  consider  a  reduced- 
order  controller  for  this  nonlinear  simulation.  Flere,  on  the  other  hand  we  determine 
the  spatial  roll-off  characteristic  of  the  controller  based  on  the  Green’s  function 
approach.  The  roll-off  characteristics  will  shed  light  on  the  spatial  resolution  of  the 
arrays  of  discrete  sensors  and  actuators  required  for  a  successful  control.  A  good 
spatial  roll-off  implies  that  relatively  few  measured  points  are  required  to  achieve  an 
effective  control  (see  Bamieh  &  Dahleh  2001).  We  refer  to  the  description  in  §2.4. 
Consider  the  same  case  in  the  numerical  simulation  for  Pr  =  7  and  Ra  =  104  and  a 
length  scale  of  the  layer  corresponding  to  ax  =  1  and  av  =  1.  Figure  6(a)  shows  the 
contour  of  the  Green’s  function  G(v,  y,  t\xp,  yp ,  tp),  which  is  the  response  temperature 
on  the  actuator  plane  z  =  0  owing  to  an  impulse  temperature  8(x—xp)8(y—yp)8(t—tp) 
on  the  sensor  plane  z*  =  0.3.  Here,  we  let  xp  =  1.5,  yp  =  1.5  and  tp  =  0,  /*  =  20 Ar  with 
At  =  10~3.  Figure  6(b)  shows  the  response  temperature  profile  as  a  function  of  x—xp 
along  the  dashed  line  designated  in  figure  6(a).  The  response  temperature  corresponds 
to  G(x,  yp,  t*\xp,  yp,  tP)  (with  t*>tp).  The  result  shows  that  the  function  has  a 
negative  minimum,  here  denoted  by  Gmin.  The  minimum  is  collocated  horizontally 
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with  the  sensor  impulse.  The  negative  temperature  generated  is  intended  to  cancel 
the  disturbance  temperature  created  by  the  impulsive  temperature.  Of  particular 
signihcance  is  how  steep  the  response  curve  (in  V-shaped)  is,  implying  that  the 
influence  zone  about  the  sensor  point  is  small.  From  figure  6(b),  the  base  width  of 
the  V-shape  curve  is  about  the  width  of  one  roll,  assuming  that  the  length  scale 
of  the  roll  does  not  differ  significantly  from  its  critical  value.  In  order  for  the  controller 
to  stabilize  the  convective  disturbance,  the  spacing  between  successive  points  in  the 
array  cannot  be  greater  than  the  effective  width  of  the  response  temperature. 

The  plots  in  figure  6(a)-6(b)  represent  a  snapshot  at  t  =  t*.  We  observed  that  as  f 
increases  from  0,  the  shape  and  width  of  the  temperature  profile  (see  figure  6b)  have 
changed  little,  but  the  magnitude  of  the  minimum  decreases  rapidly.  In  figure  6(c), 
we  show  the  change  of  the  temperature  at  the  minimum,  Gmin,  with  f .  The  large  dot 
in  figure  6(c)  denotes  the  point  corresponding  to  the  snapshot  of  figures  6(a)-6(b). 
Since  the  system  is  diffusive,  the  response  temperature  decays  monotonically  in  time, 
as  expected. 


4.  Experimental  considerations 

For  implementation  of  the  LQG  feedback  control  design,  an  experiment  of  RBC 
is  considered.  This  effort  will  be  guided  by  the  result  of  the  nonlinear  simulation, 
modified  for  air  at  room  temperature  (with  Pr  «  0.7)  as  the  working  fluid.  Although 
the  closed-loop  numerical  results  presented  earlier  in  the  paper  is  for  the  case  Pr  =  7.0 
only,  our  supplementary  analysis  completed  only  recently  at  Pr  =  0.71  has  revealed 
that  there  is  no  significant  difference  in  the  closed-loop  response  between  the  two 
Prandtl  numbers  for  the  condition  Ra  =  104. 

For  RBC,  previous  experiments  demonstrate  that  the  initial  and  onset  conditions,  as 
well  as  the  realized  convection  pattern,  are  predictable  under  controlled  experimental 
conditions  (Cross  &  Flohenberg  1993).  Complex  situations  in  applications  such  as 
variations  of  material  properties,  occurrence  of  concentration  gradient  and  solutal 
convection,  presence  of  horizontal  basic  temperature  gradient,  sidewall  effects,  defects 
in  pattern,  etc.,  are  not  included. 

In  the  experimental  apparatus,  the  upper  and  lower  walls  will  be  of  two  types  of 
material  with  a  large  range  of  heat  conductivity.  The  two  walls  have  large  aspect  ratio 
to  the  layer  depth,  and  may  have  different  thermal  boundary  conditions.  Miniature 
strain  gauge  type  heaters  will  be  strategically  placed  at  the  lower  wall  as  actuator 
(with  separation  between  heaters  determined  by  the  wavelength  of  the  pattern  to 
be  controlled).  For  air,  it  is  convenient  to  use  the  holographic  interferometry  as 
the  sensing  technique.  Such  a  sensor  can  detect  temperature  differential  to  high 
precision.  Our  LQG  controller  design  has  been  validated  using  simulated  sensor  data. 
Eventually,  for  implementation  in  the  experiment,  a  reduced-order  LQG  controller 
will  be  developed.  The  three-dimensional  pseudospectral  model  will  be  modified  to 
accommodate  the  spatial  and  temporal  dynamics  of  the  sensors  and  actuators,  guided 
by  laboratory  observations  and  the  experimental  data. 


5.  Conclusion 

The  goal  achieved  in  this  study  is  a  successful  demonstration  through  numerical 
simulations  that  a  fully  nonlinear  steady  and  preferred  state  of  convection  in  a 
horizontal  layer  of  fluid  can  be  reverted  to  the  no-motion  state  by  closed-loop 
controller  action.  The  simulated  results  here  show  the  performance  of  the  LQG 
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controller  at  Ra  =  104  and  Pr  =  7.  At  this  Ra,  the  proportional  feedback  controller  is 
ineffective  according  to  the  linear  theory.  For  even  higher  values  of  Ra,  stabilization 
is  likely  to  be  achievable  with  the  LQG  controller  by  using  higher  spatial  resolution 
in  the  simulation,  but  we  have  not  pushed  for  that  result.  The  reason  is  that  for 
realistic  modelling  at  high  Ra  the  effects  of  the  discrete  actuator  and  actuator  delay 
are  important  considerations  as  well.  Although  a  general  stability  proof  cannot  be 
inferred  from  the  nonlinear  simulation  of  a  few  initial  conditions,  the  results  do 
indicate  that  the  linear  controller  appears  quite  responsive  in  suppressing  important 
finite  disturbances. 

The  numerical  method  used  here  to  develop  the  nonlinear  plant  model  is 
pseudospectral  spatially.  The  integration  of  the  model  dynamics  equation  is  performed 
by  a  time-splitting  technique.  We  have  adopted  the  conventional  scheme  developed 
in  Marcus  (1984)  (also  see  Canuto  et  al.  1986).  However,  since  some  significant 
modification  of  the  scheme  has  been  made,  we  validate  our  fully  nonlinear  three- 
dimensional  plant  model  by  comparing  check  cases  against  published  results,  in 
particular,  from  Clever  &  Busse  (1974)  and  Busse  &  Clever  (1979).  The  agreement 
appears  reasonably  good.  Moreover,  the  direct  simulation  verifies  the  results  of  the 
weakly  nonlinear  analysis  (Or  &  Kelly  2001)  about  the  presence  of  the  controlled- 
induced  subcritical  g-hexagon  solution. 

We  have  also  examined  the  shape  function  of  the  actuator  response  by  computing 
the  Green’s  function  of  the  LQG  controller.  The  shape  of  the  actuation  temperature 
determines  the  order  of  the  horizontal  distance  between  points  of  the  sensor/actuator 
arrays  in  term  of  the  layer  gap  thickness  d.  This  information  is  of  critical  importance 
when  the  more  realistic  pointwise  sensor  and  actuator  are  used  instead  of  the 
continuous  ones. 

This  research  is  supported  by  the  United  States  Air  Force  (Grant  no.  F49620-00-1- 
0166). 
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Recent  studies  in  the  feedback  control  of  Rayleigh-Benard  convection  indicate  that  one  can  sustain  the 
no-motion  state  at  a  moderate  supercritical  Rayleigh  number  (Ra)  using  only  proportional  compensation. 
However,  stabilization  occurs  at  a  much  higher  Rayleigh  number  using  linear-quadratic-Gaussian  (LQG)  con¬ 
trol  synthesis.  The  restriction  is  that  the  convection  model  is  linear.  In  this  paper,  we  show  that  a  comparable 
degree  of  stabilization  is  achievable  for  a  fully  nonlinear  convection  state.  The  process  is  demonstrated  in  two 
stages  using  a  fully  nonlinear,  3D  Boussinseq  model,  compensated  by  a  reduced-order  LOQ  controller  and  a 
gain-schedule  table.  In  the  first  stage  a  fully-developed  convective  state  is  suppressed  through  the  control 
action  at  a  moderate  supercritical  Ra.  After  the  residual  convection  decays  to  a  sufficiently  small  amplitude,  in 
the  second  stage,  we  increase  the  Ra  by  a  large  step  and  switch  the  compensator  gains  using  the  gain-schedule 
table.  During  this  change  the  control  action  is  in  place.  Our  nonlinear  simulation  results  suggest  that  the 
nonlinear  system  can  be  stabilized  to  the  limit  predicted  by  the  linear  analysis.  The  simulation  shows  that  the 
large  Ra  jump  induces  a  large  transient  temperature  in  the  conductive  component,  which  appears  to  have  very 
small  impact  on  the  stabilization. 
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I.  INTRODUCTION 

Considerable  interest  has  been  demonstrated  in  the  sup¬ 
pression  of  convection  using  the  optimal  control  methods 
[1].  A  glimpse  at  some  important  results  can  be  found  in  [2] 
using  the  Boussinesq  model  for  Rayleigh-Benard  convection 
(RBC).  The  suppression  of  the  onset  of  RBC  by  feedback 
control  was  investigated  through  analyses,  experiments,  and 
numerical  simulations  by  a  number  of  authors.  The  results 
were  published  in  [3-11].  Potential  industrial  applications  of 
the  technique  occur  in  several  areas,  such  as  material  pro¬ 
cessing,  moulding,  and  crystal  growth.  The  outlook  is  en¬ 
couraging. 

The  proportional  feedback  control  used  by  Tang  and  Bau 
[3,7,8]  and  by  Howie  [4-6,11]  offers  a  simple,  intuitive  way 
to  control  the  system.  The  mechanism  of  suppression  corre¬ 
sponds  to  the  spatial  stabilization  of  the  unstable  perturbation 
temperature  field  by  the  temperature  control.  From  a  control 
analysis  point  of  view,  a  proportional  control  law  is  often 
considered  insufficient.  To  account  for  the  complete  spatial- 
temporal  dynamics,  linear-quadratic-Gaussian  (LQG)  syn¬ 
thesis  offers  much  better  performance,  as  demonstrated  in 
[1,9,10].  LQG  synthesis  also  provides  a  systematic  and  opti¬ 
mal  methodology  to  design  a  high-order  robust  compensator 
that  guarantees  a  certain  level  of  relative  stability  margins. 
Furthermore,  with  LQG  synthesis  there  is  an  existing  body 
of  modern  control  theory  literature  to  aid  in  the  issues  of 
controller  implementation,  systematic  order  reduction,  and 
gain  scheduling. 

We  published  two  studies  in  applying  the  LQG  synthesis 
for  controlling  RBC.  The  first  study  is  a  closed-loop  linear 
stability  analysis  [9],  in  which  detailed  stability  boundary 
curves  were  computed.  The  no-motion  state  can  be  stabilized 
up  to  14.5  times  the  critical  Rayleigh  number  Rac  for  the 
case  of  moderate  Prandtl  number  (Pr=7.0).  Beyond  this  limit 
the  LQG  controller  allows  islands  of  instability  to  form  be¬ 
low  Ra=  14.5  times  of  Rac. 


The  linear  stability  results  is  very  restricted  because  in 
real  applications  convection  has  finite  amplitudes.  In  the  sec¬ 
ond  study,  a  three-dimensional  nonlinear  plant  model  [10] 
was  developed  using  the  pseudo-spectral  method.  Using 
LQG  synthesis  of  [9],  we  demonstrated  that  a  fully  nonlinear 
initial  steady  convection  state  at  about  6  times  the  critical 
Rayleigh  number  can  be  damped  out  by  control  action.  The 
no-motion  state  is  sustained. 

In  all  the  studies  [1-10],  the  measurements  and  control 
actions  were  assumed  to  be  spatially  continuous.  No  sensor 
or  actuator  dynamics  were  included.  In  a  more  recent  study 
by  Howie  [11],  it  was  demonstrated  that  the  finite  wall  thick¬ 
ness  in  the  boundaries  adds  additional  actuator  dynamics  and 
alters  the  onset  mode  of  instability.  Here,  this  added  latency 
has  not  been  included.  This  area  is  of  interest  for  future 
investigations. 

In  this  study  our  goal  is  to  demonstrate  that  a  fully  non¬ 
linear  convection  state  can  be  suppressed,  and  the 
convection-suppressed,  no-motion  state  can  be  raised  to  a  Ra 
value  comparable  to  the  stability  limit  dictated  by  the  linear 
analysis.  A  two-stage  gain-schedule  approach  is  considered 
here.  Furthermore,  a  reduced-order  LQG  controller  has  been 
developed  for  this  task. 

II.  3D  NONLINEAR  PLANT  MODEL 

Consider  an  infinite  layer  of  Boussinesq  fluid  (see  the 
schematic  of  Fig.  1).  The  upper  wall  is  maintained  at  tem¬ 
perature  T2  and  the  lower  wall  temperature  at  T\  +  (3 AT* 
+  0C,  where  T\  >Tf  the  flag  [3  is  set  to  zero  prior  to  the  step 
jump  in  Ra  and  equal  to  one  after  the  jump;  AT*  is  the  step 
increase  in  lower-wall  temperature  corresponding  to  the  Ra 
step.  These  are  the  conductive  components  of  the  tempera¬ 
ture.  To  suppress  convection,  a  temperature  control  ff(x,y) 
is  generated  on  the  actuator  plane  assumed  to  coincide  with 
the  lower  wall.  This  is  a  perturbation  temperature. 
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FIG.  1.  Schematic  of  a  section  of  the  3D  fluid  layer. 


There  are  three  embedded  sensor  planes  in  the  layer, 
which  measure  the  interior  field  temperatures  at  three  differ¬ 
ent  levels,  z(,)  (i=  1 ,2,3).  The  perturbation  sensor  tempera¬ 
tures  from  the  conductive  component  are  denoted  as 

6*s(x,y,z(?). 

Due  to  the  step  increase  in  Ra,  a  transient  temperature  is 
generated  in  the  conductive  component  of  the  total  tempera¬ 
ture  after  the  step  increase.  The  heat  equations  for  the  static 
conductive  temperature  and  transient  conductive  tempera¬ 
ture,  and  for  the  perturbation  temperature  due  to  convection 
are  derived  from  first  principles  in  the  Appendix.  In  the  pres¬ 
ence  of  convection,  the  total  temperature  is  decomposed  into 
the  conductive  and  convective  components.  The  conductive 
temperature  has  a  transient  component  as  well.  The  nominal 
conductive  temperature  is  the  static  component.  The  conduc¬ 
tive  temperature  held  is  stable.  The  convective  temperature 
has  a  nominal  of  zero.  The  control  loop  is  designed  to  drive 
this  perturbation  component  to  zero.  In  the  nondimensional 
mathematical  form,  the  quantities  d,  d2/K,  Kid ,  K/d2, 
p(K/d)2,  and  (7^—7^)  (before  the  Ra  step)  and  (7',  +  A7^' 
-T2)  (after  the  Ra  step),  are  used  as  the  scales  of  length, 
time,  velocity,  vorticity,  pressure,  and  temperature,  where  d 
is  the  layer  thickness,  k  and  p  are  the  mean  thermal  diffusiv- 
ity  and  density  of  the  fluid. 

The  governing  nondimensional  Oberbeck-Boussinesq 
equations  are 


Pr'<?,v  =  Pr~ 

*v  X  to  +  kRa#-  V  TTf,  +  V2v, 

(1) 

II 

1 

< 

■  V  0  +  w(l  -  A,©)  +  V20, 

(2) 

dr§  =  d  20, 

(3) 

V  •  v  =  0. 

(4) 

where  \=(u,v ,w)  is  the  velocity  vector  field,  m=  V  Xv  is 
the  vorticity,  Tre=ir+v-v/2  is  the  effective  pressure  head, 
0(z, f)  is  the  transient  component  in  the  conductive  tempera¬ 
ture,  6  is  the  perturbation  temperature  due  to  convection,  and 
k  is  unit  vector  in  the  z  direction.  The  two  external  param¬ 
eters  are  the  Rayleigh  number  and  the  Prandtl  number,  re¬ 
spectively,  defined  by  R^agkT^d2’ I  vk  and  Pr=  vl  k  where 
a  is  the  coefficient  of  thermal  expansion  and  v  is  the  mean 
kinematic  viscosity.  The  continuity  equation  (4)  is  for  an 
incompressible  flow. 


The  upper  and  lower  wall  velocity  field  conditions  are 
nonpermeable  and  nonslip, 

\(x,y,0,t)  =  0,  v(x,y,  1,7)  =  0.  (5) 


The  upper  and  lower  wall  transient  conductive  temperatures 
have  a  zero  initial  condition  before  the  step  jump,  since  the 
layer  is  in  a  steady  state  heat  conduction  prior  to  the  step 
increase  of  Ra.  After  the  step  increase,  at  t=t+,  the  initial 
condition  for  ©  is  given  by 


@(z,f+)  =  - 


Art(i-z) 

(r*  +  A T*  -  t*2)  ' 


(6) 


The  boundary  conditions  are  homogeneous. 


0(0,7)  =  0,  0(1, f)  =  0.  (7) 


The  upper-wall  boundary  condition  for  the  perturbation  tem¬ 
perature  is  homogeneous.  The  lower-wall  boundary  condi¬ 
tion  for  the  perturbation  temperature  includes  the  tempera¬ 
ture  control  0c(x,y,t),  which  is  generated  by  the  actuator 
action.  The  upper  and  lower  thermal  boundary  conditions  for 
the  perturbation  temperature  are 

0{x,y,  1,7)  =  0,  0{x,y,O,t)  =  6c(x,y,t)  ■  (8) 

The  temperature  field  is  measured  on  the  three  embedded 
sensor  planes  at  the  levels  z  =  zs  (s=  1,2,3).  Since  the  con¬ 
ductive  temperature  field  in  the  layer  can  be  determined  in¬ 
dependently,  we  can  subtract  it  from  the  total  to  obtain  the 
perturbation  temperature  measurements.  The  perturbation 
sensor  temperature  measurements  are 

0s(x,y,t)  =  0{x,y,zs,t),  s  =  1,2,3.  (9) 

In  the  numerical  scheme  the  three-dimensional  (3D)  depen¬ 
dent  variables  u,  v,  w,  p,  and  0  are  expressed  by  the  follow¬ 
ing  finite  triple  sums: 


“ 
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“  ” 

u 

M/cmn 

V 

N  K  M 

V  kmn 
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(x,y,z,t)  =  Re" 

22  2 

wkmn 

7i=0  k= 0  m=—M+ 1 

p 

Pkmn 

e 

n 

ukmn 

X(z)e'^“**+mV)  >. 


(10) 


where  Re  denotes  the  real  part  of  the  sum;  ax  and  av  are  the 
fundamental  wave  numbers  in  the  x  and  y  directions,  respec¬ 
tively.  The  two  fundamental  wave  numbers  are  prescribed. 
Here,  the  number  of  total  coefficients  can  be  reduced  to  half 
by  only  including  the  coefficients  for  k  0,  because  all  the 
dependent  variables  are  real  (see  detailed  description  in 
[10]).  The  functions  Tn(z)  (h  =  0,  1 , ...)  denote  the  Cheby- 
shev  polynomials.  A  linear  coordinate  transformation  is  used 
to  convert  the  Chebyshev  function  domain  from  the  function 
domain  —  1  =£z<  1  to  our  physical  domain  0s=z<  1.  The  per- 
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turbation  temperature  control  0,  and  measurement  tempera¬ 
tures  9S  (.?=  1,2,3),  are  planar.  They  are  expanded  as  double 
sums, 


0s(z2,t) 

6c(0,t) 


1 

p 

oT' 

_ 1 

( x,y,t )  =  Re" 

K  M 

2  2 

k= 0  m=—M+ 1 

1 

cT 

s' 

_ 1 

Xe‘0caxx+mayy') 


(11) 
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0.  (z  ,  t) 
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FIG.  2.  The  control  loop  diagram. 


III.  THE  REDUCED-ORDER  LQG  COMPENSATOR 

The  LQG  synthesis  method  was  discussed  in  detail  in  [9]. 
In  the  study,  the  LQG  compensator  used  is  full  ordered,  that 
is,  the  compensator  has  the  same  order  as  the  3D  nonlinear 
convection  model.  The  controller  was  demonstrated  to  be 
effective  in  removing  the  convection  state,  reverting  the 
layer  back  to  its  no-motion  state  for  Ra=  104,  about  5.8  times 
the  critical  value.  In  this  paper,  the  full-order  compensator  is 
replaced  by  a  reduced-order  compensator.  A  methodology  for 
increasing  Ra  was  developed  without  destabilizing  the  con¬ 
trolled  state.  This  is  a  gain  scheduling  approach. 

To  design  a  reduced-order  compensator,  the  linear  Bouss- 
inesq  equations  are  expressed  in  a  state-space  form  for  each 
Fourier-decomposed  mode  [corresponding  to  wave  numbers 
(kax,may)].  The  state-space  system  includes  the  correspond¬ 
ing  measurement  (sensor)  equations  and  the  control  (actua¬ 
tor)  equations.  We  remark  that  these  are  the  linearized  equa¬ 
tions  about  the  no-motion  state,  not  about  the  nonlinear 
convective  state.  For  the  simplicity  of  notations,  the  indices 
for  the  wave  numbers  are  omitted.  It  is  clear  that  there  is  a 
distinct  set  of  state-space  equations  for  each  pair  of  wave 
numbers  for  each  Fourier  mode  for  each  Rayleigh  number. 
The  state  space  equation  is  expressed  as 

x  =  Ax  +  Bu,  (12) 

where  the  state  vector  x  consists  of  the  Chebyshev  coeffi¬ 
cients  for  the  velocity  and  the  temperature  perturbations  for 
wave  number  pair  ( kax,may );  u  is  the  Fourier  coefficient  of 
the  temperature  control  according  to  Eq.  (11).  The  sets  of 
state  space  equations  (12)  corresponding  to  different  wave- 
number  pairs  are  decoupled  in  a  linear  sense. 

From  Eqs.  (11)  and  (10),  we  denote  the  measurement  tem¬ 
perature  Fourier  coefficients  by  the  three-row  vector  z;  z 
consists  of  measurement  on  3  planes  Z\,Z2,Z 3,  which  can  be 
expressed  in  terms  of  the  state  vector  by  the  matrix  equation 

z  =  Cx.  (13) 

We  show  a  block  diagram  for  the  control  loop  between  the 
compensator  and  the  plant  model  by  a  schematic  in  Fig.  2. 
The  LQG  compensator  consists  of  a  Kalman  filter  and  an 
optimal  regulator.  The  fixed-gain  Kalman  filter  equation  and 
the  optimal  regulator  equation  are,  respectively. 


x  =  Ax  +  Bu  +  Ky(z  -  z),  m  =  -Kcx,  (14) 

where  vector  x  consists  of  the  state  estimates.  The  following 
measurement  equation  relates  the  estimates  of  measurements 
to  the  estimates  of  the  true  states, 

z  =  Cx.  (15) 

The  steady-state  Kalman  gain  vector  Ky  and  the  regulator 
gain  vector  K,  are  determined  from  two  steady-state  alge¬ 
braic  Riccati  equations  (AREs)  [1,9]. 

To  determine  the  controller  gain  Kr,  the  states  x  and  con¬ 
trol  u  are  quadratically  penalized  in  an  integral  cost  criterion. 
The  associated  weightings  are  used  in  the  controller  AREs. 
To  determine  the  filter  gain  Ky,  we  assumed  that  Gaussian 
noises  are  added  to  the  dynamics  and  the  measurements.  The 
power  spectral  densities  are  used  as  the  design  parameters  in 
the  AREs  to  form  the  gains.  The  LQG  compensator  at  each 
Ra  number  consists  of  2 (K+ 1  )M  sets  of  state  space  matrices 
(A,B,C,D)  (where  D  =  0),  and  the  same  number  of  sets  of 
gains  Ky,K(,.  There  is  no  cross  coupling  between  the  sets  of 
equations.  Each  set  is  designed  to  compensate  for  a  Fourier 
mode  corresponding  to  a  wave  number  pair.  In  the  full-order 
controller  formulation,  the  dimension  of  A  ranges  from  64  to 
128,  K  and  M  range  from  32  to  64.  The  compensator  be¬ 
comes  very  large. 

We  have  considered  two  different  approaches  to  design 
the  order-reduced  compensator.  The  first  approach  is  to  seek 
a  balance  between  the  input  and  output  relationships  of  the 
plant  model.  The  second  approach  is  to  seek  a  balance  be¬ 
tween  the  input  and  output  relationships  of  the  compensator. 
The  state  space  system  is  first  transformed  to  a  Schur  canoni¬ 
cal  form  and  then  ordered  block  by  block  in  a  balanced  re¬ 
alization  scheme  (grammian-based).  The  reduced-order  com¬ 
pensator  is  then  constructed  based  on  the  reduced-order 
dynamics  and  their  input-output  relationships. 

The  first  approach  appeared  to  be  quite  effective  when 
applied  to  controlling  the  laminar  boundary  layer  transition 
[12].  However,  when  applied  to  Rayleigh-Benard  convec¬ 
tion,  the  reduced-order  compensator  is  not  as  robust  as  it 
should  be.  Some  preliminary  simulation  results  indicate  that 
an  order  reduction  by  a  factor  of  2  is  possible.  Further  reduc¬ 
tion  makes  the  closed-loop  response  unstable. 


046302-3 


A.  C.  OR  AND  J.  L.  SPEYER 


PHYSICAL  REVIEW  E  71,  046302  (2005) 


The  first  approach  is  not  robust.  Our  ultimate  goal  is  to 
use  a  reduced-order  compensator  to  control  a  full-order  non¬ 
linear  plant  model.  The  reason  that  we  perform  a  balance 
realization  on  the  linearized  plant  model,  and  then  truncate 
the  plant  model,  is  because  we  want  to  use  the  reduced-order 
linear  plant  model  to  design  the  reduced-order  compensator. 
In  the  second  approach,  we  do  not  truncate  the  plant  model. 
We  first  design  a  full-order  compensator  for  the  full-order 
plant  model,  just  as  in  [9].  Then,  we  seek  to  perform  a  bal¬ 
ance  realization  on  the  compensator  instead,  and  then  trun¬ 
cate  compensator  model.  By  doing  that  we  have  to  compute 
the  full-order  Kalman  gains  K,  and  full-order  regulator  gain 
Ke  based  on  the  full-order  plant  model  with  state  space  (A, 
B,  and  C).  Then,  we  perform  a  balance  realization  on  the 
compensator  input  and  output.  The  compensator  inputs  the 
measurement  vector  z  and  outputs  the  control  vector  u. 

From  Eq.  (14),  we  obtain  the  following  state  space  equa¬ 
tions  for  the  full-order  compensator. 


x  =  Acx  +  Brz, 
u  =  Cfx  +  Dt,z . 


(16) 


The  subscript  “c”  denotes  the  compensator,  where  the  state 
space  matrices  are  given  by 


Ac  =  A-  BKC  -  K7C  +  KyDKc, 

(17) 

B(=K/;  Cc  =  -  Kc,  D£  =  0, 

where  I),  =0.  To  reiterate,  in  the  first  approach  the  balance 
realization  is  performed  on  (A,B,C).  In  the  second  ap¬ 
proach,  the  balance  realization  is  performed  on 
(A(„Bt„Cc,D  c)  instead.  It  is  worth  noting  that  the  dynamics 
of  Ac  can  be  distinctively  different  from  those  of  A.  The 
compensator  dynamics  now  depends  on  the  gain  vectors  as 
the  input  and  output  matrices. 

A  grammian-based  balance  realization  is  just  one  way  to 
perform  the  balancing.  A  typical  grammian-based  routine  is 
by  solving  the  Lyapunov-type  equations  using  the  polyno¬ 
mial  methods.  The  polynomial  routines  are  not  effective  for 
large  matrices.  Here,  we  consider  an  approximate  method. 
First,  we  apply  a  similarity  transformation  x  =  Px  h  in  which 
the  transformation  P  diagonalizes  Ac.  This  transformation 
converts  (Ac,Bc,Cc,Dc)  into  (Acl,Bcl,Ccl,Dcl).  If  there  are 
nondistinct  eigenmodes,  then  we  have  to  use  the  Schur  ca¬ 
nonical  transformation  instead.  The  transformation  is 


Acl=PAfP,  Bcl  =  P-1BC, 
Cci  =  CcP,  Dcl  =  Dc, 


(18) 


where  Acl  is  now  in  a  diagonal  form.  We  denote  the  diagonal 
entries  of  this  matrix  by  A,-  (i=l,2, ...  ,N).  We  allow  the 
diagonal  entries  to  be  expressed  in  complex  conjugate  pairs. 

Second,  we  define  a  score  parameter  for  each  compensa¬ 
tor  mode,  denoted  as  cr,  (i=  1,2,... ,  A'),  where 

ai =  |b,| |c,j/| A,|,  i  =  1,2,  ...  ,N,  (19) 

b,  is  the  zth  row  of  Bel  and  c,  is  the  zth  column  of  Ccl.  Then 
we  perform  a  reordering,  in  which  the  eigenmodes  in 


(Afl,Bcl,Ccl,Dcl)  are  ranked  according  to  the  magnitude  of 
the  score  parameter.  The  modes  are  arranged  in  the  order  of 
importance.  The  unimportant  modes  at  the  end  will  be  trun¬ 
cated.  A  truncated,  reduced-order  set  is  denoted  by 
(Ac2r,BC2r,CC2r,Dc2r),  with  the  state  vector  xlr  (subscript 
“r”  stands  for  reduced-order  set). 

The  reduced-order  compensator  is  described  by 


^lr  ^r.2A  !  r  ^ 

u  =  Cc2rxlr+  Dc2rz. 


(20) 


IV.  RESULTS 

The  preferred  convection  pattern  is  steady,  two- 
dimensional  (2D)  finite-amplitude  rolls.  We  start  at  Ra=104 
and  use  the  steady  convection  rolls  as  the  initial  condition. 
Although  the  preferred  wavenumber  of  the  rolls  is  about  3.1, 
in  our  simulation,  the  fundamental  wavenumbers  are  chosen 
to  be  at  a:=  1.0  and  ay=  1.0.  So,  the  steady  rolls  occur  as  the 
longitudinal  rolls  (parallel  to  y  axis)  in  the  second  Fourier 
harmonic.  Three-dimensional  zero-mean  Gaussian  noises 
(zero  mean,  standard  deviation  equal  to  10-6)  are  added  to 
each  coefficient  in  the  initial  fields.  This  ensures  that  if  3D 
instabilities  occur,  there  exists  perturbations  for  them  to 
grow.  The  full-order  nonlinear  plant  model  has  a  resolution 
of  2  X  128  Chebyshev  modes  in  the  vertical  direction  (for 
both  vertical  velocity  and  temperature)  and  32  X  32  Fourier 
modes  in  the  horizontal  direction.  Note  that  the  horizontal 
velocity  is  obtainable  from  the  vertical  velocity  through  the 
continuity  equation.  In  the  reduced-order  compensator,  only 
eight  Chebyshev  modes  (complex)  are  retained.  The  same 
number  of  Fourier  modes  is  used  in  the  horizontal  direction. 
The  analysis  shows  that  8  Chebyshev  modes  is  the  minimum 
for  effective  control. 

A.  Closed-loop  simulation  prior  to  the  Ra  increase 

The  first-stage  of  simulation  is  performed  using  the  first 
set  of  gains  in  the  gain-schedule  table.  The  simulation  is 
performed  until  the  residual  convection  is  sufficiently  small. 
The  simulation  shows  that  a  time  period  of  0.4  unit  is  suffi¬ 
cient  for  the  purpose.  With  a  time  step  of  Ar= 0.0004,  we 
sample  the  output  every  five  time  steps.  The  maximum  tem¬ 
perature  control  recorded  occurs  right  at  the  start,  i.e.,  at  1 
=  0.004.  The  temperature  control  produces  a  sharp  peak  in 
the  transient  phase  of  roughly  0.15  time  unit  wide.  This  peak 
has  a  maximum  roughly  of  1.3  times  the  conductive  tem¬ 
perature  difference  imposed  across  the  upper  and  lower 
walls,  that  is,  A T.  This  1.3AT  may  be  considered  high  for 
experimental  implementations.  It  all  depends  on  the  satura¬ 
tion  limit  of  the  actuator.  If  saturation  occurs,  then  we  have 
to  reduce  the  Ra  value  of  the  initial  condition. 

The  question  of  using  a  linear  controller  to  control  a  non¬ 
linear  system  needs  an  explanation.  In  fact,  the  linear  system 
consists  of  only  the  one  Fourier  mode.  Therefore,  in  the  lin¬ 
ear  sense  the  controller  is  a  single-mode  controller  [3-6,9]. 
The  nonlinear  system  consists  of  an  infinite  number  of  Fou¬ 
rier  modes,  all  except  one  are  generated  by  the  nonlinear 
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FIG.  3.  Closed-loop  responses 
on  three  test  cases. 


terms.  Other  than  the  Fourier  modes,  there  always  exists  a 
mean  temperature  field  as  well.  The  linear  controller  controls 
each  Fourier  mode,  as  well  as  the  mean  field,  individually. 
Even  though  each  mode  is  treated  as  linear  by  the  controller, 
a  linear  controller  has  been  proven  effective.  In  fact,  it 
should  be  effective,  as  long  as  the  higher-order  terms  in  the 
amplitude  expansion  corresponding  to  a  particular  Fourier 
mode  is  either  stabilizing  or  insignificant.  Therefore,  the  tem¬ 
perature  control  containing  a  horizontal  mean  component  is 
obvious  but  was  not  explicitly  mentioned  in  [9].  This  mean 
component  is  reasonably  large  too.  At  the  peak  of  the  tem¬ 
perature  control,  this  component  is  about  -0.35A T.  It  tends 
to  zero  as  convection  is  suppressed. 

The  Nusselt  number  Nu  is  defined  as  the  ratio  of  the  total 
(convective  and  conductive)  to  the  purely  conductive  heat 
transfer.  This  parameter  is  a  measure  of  the  degree  of  con¬ 
vective  activities.  In  the  closed-loop  response,  it  is  also  a 
measure  of  the  effectiveness  of  the  controller.  In  Figs.  3  and 
4,  the  Nusselt  numbers  at  the  lower  wall  (Nu1;  solid  line)  and 
upper  wall  (Nu2,  dashed  line)  are  plotted  as  functions  of 
time.  When  Nu1  =  Nu2=1.0,  the  layer  is  purely  conductive, 
i.e.,  there  is  no  convection  motions.  However,  the  conductive 
state  can  be  transient.  When  Nu,  >Nu2,  relative  to  conduc¬ 
tion,  net  heat  is  pumped  into  the  layer  via  the  lower  wall. 
This  situation  is  observed  at  the  beginning  of  the  closed-loop 
simulation  (for  time  less  than  0.07)  in  every  case,  indicating 
that  the  controller  starts  by  pumping  net  heat  into  the  layer. 
When  Nu2>Nu1;  net  heat  (relative  to  conduction)  is  ex¬ 
tracted  from  the  layer.  It  is  possible  for  the  transience  to 
cause  Nui  <  1.0  as  convection  is  damping  out,  as  seen  in  the 
conductive  case  (see  Fig.  4).  At  small  convection  amplitude 
(as  the  Nusselt  number  approaches  1.0),  the  nominal  simula¬ 
tion  (see  Fig.  4)  shows  that  the  upper-wall  Nusselt  number 


has  value  greater  than  1 .0  and  the  lower-wall  Nusselt  number 
has  value  less  than  1.0.  This  linear  behavior  shows  that  the 
controller  action  is  to  extract  net  heat  (relative  to  conduction) 
from  the  layer,  thereby  reducing  the  convection  energy. 

B.  Gain  scheduling,  Ra  stepping  and  the  transient  response 
of  the  conductive  temperature 

The  gain  scheduling  is  a  means  to  allow  the  filter  and 
control  gains  to  vary  during  the  closed-loop  simulation.  Typi¬ 
cally,  a  gain  table  is  constructed.  The  Ra  range  is  divided 
into  intervals.  Each  interval  corresponds  to  a  new  set  of 
gains.  In  the  simulation,  as  Ra  is  increased  through  the  inter¬ 
vals  in  time,  the  gains  are  switched  to  preserve  the  control 
performance.  For  implementation,  it  is  more  convenient  to 
increase  Ra  by  steps.  A  step  jump  causes  transience  in  the 
conductive  temperature  component.  This  may  in  turn  cause 
perturbations  in  the  convective  field.  In  this  study,  we  con¬ 
sider  only  a  two-interval  gain-table,  since  each  set  of  gains  is 
extremely  large  even  for  the  reduced-order  controller  [32 
X32  sets  of  reduced-order  state  space  (A,B,C,D)].  We  use 
a  large  step  increase  in  Ra  to  assess  the  performance.  The 
key,  however,  is  not  to  destabilize  the  closed-loop  system. 

In  the  simulation  we  consider  a  step  increase  of  Ra  from 
104  to  2X104.  The  simulation  following  the  step  at  r=t+ 
=  0.4  unit  continues  for  0.4  time  unit.  The  total  simulation 
time  is  0.8  unit. 

Before  performing  the  nominal  simulation,  we  consider 
three  test  cases  (shown  in  the  three  panels  of  Fig.  3).  The 
purpose  of  the  tests  is  for  sanity  checks  and  also  for  a  better 
understanding  of  the  closed-loop  response  in  general.  In  the 
test  case  1  (see  Fig.  3,  upper  panel),  the  control  action  is 
disconnected  right  after  the  Ra  step  increase  at  r=r+=0.4. 
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TIME  IN  DIFFUSIVE  UNIT 


FIG.  4.  Closed-loop  upper  and 
lower  Nusselt  numbers. 


The  remaining  simulation  has  no  feedback  control.  The  key 
observation  is  that  it  takes  a  long  duration  for  convection  to 
build  up  again.  Eventually,  the  layer  has  reverted  the  sup¬ 
pressed  state  to  full-blown,  finite-amplitude  steady-state  of 
convection.  The  behavior  is  indicated  by  the  Nu2  (dashed 
line)  and  Nu ,  (solid  line).  The  curves  converge  to  its  steady 
state  value  at  about  3.0. 

In  the  test  case  2  we  maintain  a  closed  loop  throughout 
the  simulation.  However,  we  do  not  switch  the  gains  at  the 
step  increase  of  Ra.  Figure  3  (the  middle  panel)  shows  that 
the  original  set  of  gains  is  ineffective  to  maintain  the  no¬ 
motion  state.  Even  worse  than  the  test  case  1,  the  closed  loop 
response  displays  finite-time  divergence.  In  this  case,  the 
system  does  not  settle  to  a  steady  state.  It  blows  up. 

A  set  of  gains  designed  at  a  higher  Ra  is  not  necessarily 
effective  to  control  the  layer  at  a  lower  Ra.  In  the  test  case  3, 
we  use  the  second  set  of  gains  (designed  at  Ra=2  X  104)  to 
control  convection  before  the  Ra  step  increase  (the  steady 
state  at  Ra=  104).  The  result  (Fig.  3,  lower  panel)  shows  that 
the  closed  loop  response  is  divergent. 

Finally,  we  turn  to  the  nominal  case  of  the  two  Ra-step 
simulation.  In  this  case,  we  use  the  gains  according  to  the 
gain  table,  show  the  closed-Loop  response  in  Fig.  4.  At  t 
=  f+=0.4,  we  see  a  very  small  kink  in  the  solid  curve  corre¬ 
sponding  to  Nu,,  due  to  the  transient  conductive  temperature. 
The  important  observation  is  that  the  transient  conductive 
temperature  has  very  little  impact  on  the  residual  convection. 

Unfortunately,  the  Nusselt  number  plots  do  not  show  the 
detailed  transience  of  the  residual  convection  on  the  scale 
used  to  plot  Fig.  4.  It  is  of  interest  to  show  the  control  tem¬ 
perature  6C  in  the  closed-loop  response.  Referring  to  Eq. 
(11),  we  define  the  plotted  quantity  ||0f||  as 


(KM  1/2 

2  2  WMWo,*)f  ■  (2i) 

k= 0  m=-M+ 1  J 

This  measure  is  plotted  in  a  semilogarithm  scale  (see  Fig.  5). 
The  regions  of  steep  slopes  indicate  locations  where  a  cross 
over  zero  has  occurred.  The  maximum  of  the  mean  tempera¬ 
ture  control  occurs  at  the  beginning  of  the  simulation.  The 
maximum  value  is  1.27.  Near  r=0.4,  a  distinct  jump  is  ob¬ 
served.  This  is  due  to  the  conductive  temperature  transience 
caused  by  the  step  increase  in  Ra.  The  dominant  coefficient 
(greater  than  96%  of  the  norm,  corresponds  to  the  coefficient 
for  k= 3  and  m  =  0,  as  expected.  This  coefficient  represents 
the  longitudinal  convection  rolls  at  a  wave  number  3.0,  as 
prescribed.  All  the  coefficients  for  m  A  0  vanish.  The  results 
suggest  that  a  large  step  increase  of  Ra  is  benign.  The  large 
transient  conductive  temperature  does  not  destabilize  the 
closed-loop  system. 

In  Fig.  6,  we  plot  the  transient  conductive  temperature, 
©(z,/j  =  ©(z,f)-©„(z)  at  several  time  points,  at  r=0.4,  0.48, 
0.56,  and  0.72.  The  quantity  converges  to  zero  fairly  rapidly. 

V.  CONCLUSIONS 

In  this  study  two  important  steps  have  been  advanced  us¬ 
ing  a  robust  control  synthesis  for  the  problem  of  suppression 
of  convection:  (i)  a  reduced-order  compensator  and  (ii)  a 
simulation  based  on  a  two-step  gain-schedule  table.  For  com¬ 
pensator  order  reduction,  the  simulation  results  show  that 
eight  Chebyshev  modes  is  the  minimum  for  an  effective  per¬ 
formance.  In  stepping  the  Ra  by  a  large  increment  using  the 
gain  schedule  table,  an  interesting  observation  is  that  a  large 
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FIG.  5.  Norm  of  closed-loop 
control  temperature  coefficients. 


Ra  step  is  accompanied  by  a  large  transient  conductive  tem¬ 
perature,  but  does  not  cause  a  significant  transient  response 
in  the  convective  fields.  As  long  as  the  residual  convection 
prior  to  the  step  jump  of  Ra  remains  small,  the  destabilizing 
effects  of  the  transient  conductive  temperature  is  insignifi¬ 
cant.  In  this  study  the  two-step  Ra  gain  table  used  for  dem¬ 
onstration  purpose.  A  more  refined,  multistep  gain  table  may 
be  necessary  for  laboratory  implementations.  In  conclusion. 


to  suppress  a  fully  nonlinear  convection,  the  gain-schedule 
approach  is  a  means  to  achieve  a  performance  up  to  the 
linear  stability  limit. 

ACKNOWLEDGMENT 

This  research  is  supported  by  the  United  States  Air  Force 
(Grant  No.  F49620-00-1-0166). 


0.2 


0.1 


0  ** ;  - 


s  -0-1 

w 

® 

I 

©  -0.2  - 

-0.3  - 

-0.4 


•  r 

:l  / 

:i  • 

V 

N. 

% 

> 

** 

‘  * 

*»» 

- "  ' 

-  t=0.40 

-  -  t=0.48 
t=0.56 
t=0.72 

-0.5 


0.1  0.2  0.3  0.4 


0.6  0.7  0.8  0.9 


FIG.  6.  Time  plot  of  the  tran¬ 
sient  conductive  temperature. 


_ i _ 

0.5 

Z-AXIS 


046302-7 


A.  C.  OR  AND  J.  L.  SPEYER 


PHYSICAL  REVIEW  E  71,  046302  (2005) 


APPENDIX:  DERIVATION  OF  THE  MODEL  EQUATIONS 

The  dimensional  form  of  the  Boussinesq  system  consists 
of  the  momentum  equation, 


0*(z,i)  =  0»  +  0(u), 


©;s(z)  =  (1*1  +  Ar*)  -  (r*  +  AT*  -  T*2)z. 


(A6) 


9,y'  +  v*  •  Vv*=-V  ( p  +  pgz)  -  ag(T*0  -  7')k  +  vV2y* , 

P 

(Al) 

where  in  the  buoyant  force  term,  T0  is  a  reference  tempera¬ 
ture  (at  z=0).  The  heat  equation  is 

9,T*  +  y*-  VT*=kV2T".  (A2) 

The  velocity  field  y* ={u*  ,v*  ,w*)  satisfies  the  incompress¬ 
ibility  condition 

V  •  v*  =  0.  (A3) 

We  use  (x,y,z)  and  t  to  denote  dimensional  as  well  as  non- 
dimensional  Cartesian  coordinates  and  time  for  simplicity. 
The  asterisk  highlights  dimensional  variables.  For  nondimen- 
sionalization  the  layer  thickness  d  and  d2l  k  are  used  for 
scaling  length  and  time;  v  and  k  denote  the  viscous  and 
thermal  diffusivity,  respectively;  p  is  the  mean  density  of  the 
fluid;  a  is  the  coefficient  of  thermal  expansion  and  g  is  the 
gravitational  acceleration.  The  upper  (at  z=l)  and  lower  (at 
Z=0)  walls  of  the  layer  are  prescribed,  respectively,  at  tem¬ 
peratures  T*2  and  T\  with  T*  >  T*2. 

The  fluid  temperature,  T*(x,y  ,z,t),  is  decomposed  into  a 
conductive  temperature  ©  *{z,t)  and  a  convective  (perturba¬ 
tion)  temperature  9{x,y,z,t), 

T*(x,y,z,t)  =  0*(z,f)  +  9*{x,y,z,t).  (A4) 

For  most  studies  in  RBC,  a  static  (nominal)  state  is  used, 

<d*  =  ®*ss=T\-{T\-T*2)z.  (A5) 


The  transient  conductive  temperature  0(z, f)  is  governed  by 

<?,©*  =  *4©*,  (A7) 

subjected  to  initial  and  boundary  conditions 

©*(z,0)  =  -  AT*(1  -  z),  0*(O,f)  =  0,  0*(U)  =  O. 

(A8) 

Next,  we  incorporate  the  decomposed  temperatures  into  the 
governing  equations  (Al)  and  (A2).  The  conductive  tempera¬ 
ture  is  balanced  by  pressure.  Equation  (Al)  gives 

9,y*  +  v*  •  V  v*  =  -  -  V  (tt  )  +  ag  0  k  +  vV2v*, 

P 

(A9) 

T T  =p  + pgz  +  I  [0*(O,f)  -  0*(z',f)]Jz'. 

J  o 

This  equation  has  the  same  form  regardless  of  the  conductive 
temperature.  Separating  the  total  temperature  into  the  con¬ 
ductive  and  perturbation  components,  the  heat  equation  be¬ 
comes 

9,6"  +  \  ■  V  6*  =  -  w*9z 0*  +  kW26* .  (A10) 

Equations  (A3),  (A7),  (A9),  and  (A10)  provide  the  governing 
equations. 

Introducing  the  nondimensional  scaling,  the  full  nondi- 
mensional  governing  equations  of  the  plant  model  are 

Pr“1((?fv  +  v  •  Vv)=-  V-7r+Ra0k  +  V2v,  (All) 
V  •  v  =  0,  (A12) 

9,6+  v  •  V  6  =  w{  1  -  <7,0)  +  V20,  (A13) 


For  Ra  >  Rac,  the  steady  perturbation  temperature  0  V  0 
does  not  vanish. 

Here  the  full  set  of  equations  are  rederived  for  the  case 
when  the  conductive  temperature  is  not  in  a  steady  state.  The 
gain  schedule  algorithm  is  applied,  when  the  Rayleigh  num¬ 
ber  is  increased  incrementally  in  steps.  Consider  the  conduc¬ 
tive  state  to  be  static  as  described  by  Eq.  (A5)  at  t=t~.  In  the 
closed-loop  response  0(x,y,z,t)~ 0  as  convection  is  sup¬ 
pressed.  At  t=t+  the  lower  wall  temperature  is  increased  to 
Tl  +  AT*,  where  r~t+.  The  conductive  temperature  after  the 
jump  is  given  by 


<?r@  =  4@. 


(A14) 


In  the  gain  schedule  approach,  the  step  increase  in  Ra  im¬ 
parts  a  discontinuity  in  the  thermal  field  of  the  system.  This 
condition  translates  an  initial  condition  for  the  conductive 
temperature  at  t=t+. 


©(z,n=- 


AT*(1  -z) 

(T*  +  AT*  -  T*2) ' 


(A15) 


Thus  a  transient  response  follows  each  jump  ©  where  0 
— *  0  as  the  conductive  temperature  approaches  steady  state. 
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The  robustness  of  control  is  a  requirement  to  maintain  a  fluid  layer  at  conductive  equilibrium  heated  to  a 
highly  supercritical  condition.  Robustness  determines  how  much  uncertainties,  or  design  parameter  mis¬ 
matches,  can  be  tolerated.  Both  linear  stability  analysis  and  three-dimensional  fully  nonlinear  simulations  are 
used  for  the  study  of  the  linear  quadratic  Gaussian  (LQG)  controller.  The  parameter  mismatches  from  the 
nominal  conditions  are  introduced  into  the  plant  model,  while  the  LQG  compensator  assumes  nominal  condi¬ 
tions.  The  mismatches  arise  from  boundary  properties,  actuator  lag,  sensor  level  uncertainty,  and  wall  thick¬ 
ness,  as  well  as  from  the  major  parameters  such  as  Prandtl  number,  Rayleigh  number,  wave  number,  and 
truncation  number  in  the  reduced-order  model.  The  results  suggest  that  the  LQG  compensator  action  can 
preserve  closed-loop  stability  at  over  ten  times  the  critical  Rayleigh  number,  provided  that  the  mismatches  in 
the  sensor  level  and  wall  thickness  are  small.  Mismatches  in  the  Prandtl  number  and  wall  material  properties 
have  little  impact.  Mismatches  in  Rayleigh  number  and  wave  number  are  relatively  benign  compared  with  the 
sensor  and  thickness  parameters.  Techniques  for  measuring  the  plant  output  temperature  at  multiple  levels  with 
sufficient  accuracy  may  be  an  implementation  challenge. 

DOI:  10.1 103/PhysRevE. 73. 046307  PACS  number(s):  47.20.Bp 


I.  INTRODUCTION 

Fluid  flow  control  has  become  a  rapidly  growing  area  of 
research,  that  combines  traditional  knowledge  of  fluid  dy¬ 
namics  and  numerical  computations  and  modern  estimation 
and  control  theory.  A  particular  well-known  problem  that  has 
been  pursued  by  investigators  with  interest  in  recent  years  is 
the  convection  suppression  in  a  layer  of  fluid  at  supercritical 
condition  (a  simple  two  degree-of-freedom  analogy  is  con¬ 
trolling  an  inverted  pendulum  about  its  vertical  position). 

The  goal  of  analysis  is  to  develop  a  robust  compensator 
design  for  convection  suppression  using  a  modern  control 
synthesis  approach.  The  Boussinesq  model  of  Rayleigh- 
Benard  convection  has  been  investigated  quite  vigorously 
because  the  experimental  configuration  can  be  modeled  and 
simulated  with  high  fidelity.  In  Refs.  [1,2]  a  purely  propor¬ 
tional  feedback  control  law  is  shown  to  produce  encouraging 
results,  even  though  the  control  law  is  well  known  for  its 
limitations.  Another  limitation  is  the  choice  and  placement  of 
temperature  sensors,  Howie  [1,3]  measured  the  mean  tem¬ 
perature,  averaged  over  the  fluid  depth  (by  shadowgraph 
technique);  Tang  and  Bau  [2,4]  measured  the  temperature 
only  at  a  single  level  along  the  fluid  depth.  Neither  has  pro¬ 
vided  temperature  measurements  sufficient  to  characterize 
the  disturbance  mode  shapes  along  the  direction  of  the  layer 
depth.  The  mixed  findings,  discrepancies  between  theoretical 
and  experimental  results  reported  by  the  authors,  may  very 
well  be  due  to  the  lack  of  robustness  of  control,  the  simple 
feedback  control  law,  as  well  as  the  crude  sensor  outputs  and 
lack  of  observability. 

Measuring  disturbance  temperatures  at  multiple  levels  im¬ 
proves  the  system  observability  significantly.  In  our  model, 
knowledge  of  temperatures  at  three  interior  levels  of  the  fluid 
layer  is  assumed.  This  is  a  significant  assumption.  For  an 
actual  laboratory  implementation,  the  remote  infrared  (IR) 


sensing  technique  is  probably  the  only  viable  way  to  acquire 
such  knowledge.  Furthermore,  for  a  complete  horizontal  cov¬ 
erage  of  the  fluid  layer,  it  will  probably  require  a  high- 
frequency  scanning  and  sensing  technique.  This  technique  is 
available  for  many  applications  but  is  not  investigated  in  this 
study. 

To  improve  the  robustness  of  control  we  have  developed  a 
linear  quadratic  Gaussian  (LQG)  control  synthesis  compen¬ 
sator  [5].  The  LQG  design  with  loop  transfer  recovery  (LTR) 
is  well  known  and  proven  to  give  robust  performance  in 
many  applications.  The  LQG-LTR  compensator  is  applied  to 
a  three-dimensional,  fully  nonlinear  model  [6]  developed 
based  on  a  Boussinesq  system  of  equations  treated  by  a  spec¬ 
tral,  time-splitting  technique  (see  Marcus  [7]).  The  results 
show  that  the  linear  compensator  is  effective  in  damping  out 
an  initial  state  of  finite-amplitude  convection  at  a  higher  Ra 
than  that  reported  in  Refs.  [1,4].  The  linear  compensator  is 
modally  distributed.  Although  the  compensator  built  for  each 
Fourier  mode  is  linear,  controlling  all  Fourier  harmonic 
modes  is  nonlinear  in  the  sense  that  all  these  modes  are  gen¬ 
erated  by  the  advective  nonlinearities.  Two  further  significant 
advances  in  the  compensator  design  are  reported  in  Ref.  [8], 
These  are  (i)  the  successful  building  of  an  order-reduced 
LQG  compensator,  and  (ii)  the  developing  of  a  gain  schedule 
algorithm. 

This  paper  aims  at  the  robustness  of  the  LQG  compensa¬ 
tor.  To  enhance  the  fidelity  of  the  existing  plant  models  (lin¬ 
ear  model  for  stability  analysis  and  fully  nonlinear  and  three- 
dimensional  for  simulations),  we  incorporate  finite  walls  to 
replace  the  isothermal  boundary  condition.  A  finite-wall 
model  is  studied  by  Howie  [9].  In  the  parameter  region  stud¬ 
ied  the  closed-loop  behavior  has  an  interesting  oscillatory 
mode.  As  for  our  LQG  compensator  model,  both  the  ideal¬ 
ized  and  finite-wall  versions  are  used.  To  avoid  confusion, 
we  denote  the  compensator  with  idealized  boundary  condi- 
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sensor  plane  j  actuator  plane 


FIG.  1.  The  sensor  and  actuator  location  in  the  fluid  layer. 

tions  [5]  as  Cl  and  the  compensator  with  finite  conductive 
walls  as  C3.  Also,  in  addition  to  using  one  actuator  plane  at 
the  bottom  wall  (referred  as  a  single-plane  actuator),  an  ad¬ 
ditional  actuator  plane  at  the  upper  wall  (referred  as  a 
double-plane  actuator)  is  used.  For  the  evaluation  of  control¬ 
ler  robustness,  in  our  previous  study  [5],  the  plant,  like  the 
compensator,  is  set  at  the  nominal  parameters. 

In  Ref.  [5],  the  Nyquist  criterion  provides  the  indicator  for 
robustness.  The  Nyquist  criterion  gives  the  total  margin. 
Moreover,  no  detailed  investigation  of  the  stability  margin  as 
a  function  of  the  wave  number  is  provided.  In  this  study, 
effects  that  impact  robustness  will  be  sorted  out  individually, 
as  a  function  of  wave  number.  The  stability  margins  are  ob¬ 
tained  with  respect  to  the  Rayleigh  number,  actuator  lag, 
sensor-plane  depth  uncertainty,  wall  thickness  uncertainty, 
and  single-plane  versus  double-plane  actuators.  Besides 
computing  the  stability  limiting  curves  based  on  linear  mod¬ 
els,  we  also  introduce  a  singular-value  bound  as  a  more  con¬ 
servative  criterion  for  robustness,  which  is  applicable  to 
multi-input-multi-output  (MIMO)  plant  systems.  Besides  the 
linear  stability,  we  also  investigate  robustness  using  a  fully 
nonlinear,  three-dimensional  (3D)  simulation  tool. 

II.  MATHEMATICAL  FORMULATION 
A.  Plant  model 

In  Fig.  1  we  show  a  schematic  of  a  section  of  an  infinite 
layer  of  fluid,  bounded  by  two  finite-conductive  walls  with 
finite  thickness.  The  fluid  layer  has  thickness  d*  [in  this  pa¬ 
per  an  asterisk  denotes  dimensional  variables  or  compensator 
parameters  (at  nominal  conditions)].  The  outer  surfaces  of 
the  lower  and  upper  walls  are  prescribed  at  temperatures  T2 
and  7] ,  respectively.  The  upper  and  lower  wall  of  thicknesses 
are  du  and  d(.  For  nondimensional  scalings  we  use  the  exter¬ 
nal  temperature  difference  AT*  =  T2-T\  as  the  temperature 
scale  and  the  fluid  layer  thickness  <F  as  the  length  scale.  In 
addition,  the  fluid  thermal  diffusive  time  scale  d*2/  k  is  used 
for  scaling  time,  where  k  is  the  fluid’s  thermal  diffusivity. 
The  equilibrium  conductive  temperature  modified  by  the 
presence  of  the  finite  walls  is  derived  in  the  Appendix.  In  the 
nondimensional  form,  the  fluid  is  governed  by  the  Oberbeck- 
Boussinesq  (OB)  equations, 

Pr“ 1  9,\  =  Pr“ 1  v  X  to  +  kRa0-  V7re  +  V2v,  (1) 
9,0=  -  v  •  V0+  w  +  V20,  (2) 


V  •  v  =  0,  (3) 

where  \  =  (u,v  ,w)  is  the  velocity  vector  held,  u>= V  X  v  is  the 
vorticity,  ire=  ir+\- v/2  is  the  effective  pressure  head,  0  is 
the  perturbation  temperature  with  respect  to  the  equilibrium 
conductive  state.  The  z  direction  is  directed  upward  and  k  is 
the  unit  vector. 

The  plant  system  has  two  external  nondimensional  param¬ 
eters:  the  Rayleigh  number  and  the  Prandtl  number.  The 
Prandtl  number  is  Pr  =v  / k*,  where  v  is  the  kinematic  vis¬ 
cosity.  The  Rayleigh  number  is  defined  in  terms  of  the  outer- 
wall  temperature  difference  and  the  fluid  layer  thickness, 
Ra=a*gAT*(d*)3/v*K*,  where  a  is  the  coefficient  of  ther¬ 
mal  expansion.  All  the  material  properties  are  assumed  con¬ 
stant.  Note  that  Ra  here  is  defined  based  on  the  outer  wall- 
to-wall  temperatures  AT*  but  fluid  thickness  d*.  The 
effective  Ra^  should  be  defined  based  on  the  fluid  boundary 
temperatures.  In  the  Appendix,  Ray  is  derived  from  Ra,  by 
replacing  dt  with  d  h,  where  h  is  a  nondimensional  factor. 
Equation  (3)  is  the  continuity  equation  for  incompressible 
flow. 

In  addition  to  the  OB  equations  describing  the  fluid  layer, 
the  lower  and  upper  wall  walls  are  governed  only  by  the 
heat-transfer  equations, 


9, 0t  =  k(V20(, 

(4) 

9,0u=kuV2Ou, 

(5) 

where  Kf  =  K*(/K  and  ku=kJk. 

The  kinematic  boundary  conditions  applied  to  the  velocity 
held  in  the  fluid  layer,  which  are  defined  as  the  nonperme- 
able  and  nonslip  conditions, 

v(x,y,0,f)  =  0,  v(x,y,  l,r)  =  0.  (6) 

The  fluid-wall  interfacial  thermal  boundary  conditions  are 
dynamic,  not  prescribed.  They  are  matched  conditions  of 
temperatures  and  heat  flux, 

0(x,y,O,f)  =  0t(x,y,0,t), 

9z0(x,y,z,t)  Lo  =  k(d.0f(x,y,z,t)  Uo  (7) 

and 

0(x,y,  l,l)  =  9u(x,y,\,t), 

9z0{x,y,z,t)  |-=i  =  kudz0u(x,y,z,t)\z=  i  ■  (8) 

The  parameters  kf  and  ku  are  the  thermal  conductivity  ratios 
of  the  upper  and  lower  wall  materials  defined  in  the  Appen¬ 
dix,  respectively. 

The  thermal  actuators  correspond  to  the  prescribed  outer 
wall  thermal  boundary  conditions, 

0f{x,y,- de,t )  =  u((x,y,t),  (9) 

where  Uf(x,y,t)  is  the  temperature  control  input  imposed  at 
lower  wall,  and 


046307-2 


ROBUST  CONTROL  FOR  CONVECTION  SUPPRESSION... 

9u(x,y,  1  +  du,t)  =  uu(x,y,t) ,  (10) 

where  uu(x,y,t)  is  the  temperature  control  input  imposed  at 
upper  wall  in  case  of  two  actuator  planes.  The  control  inputs 
to  the  plant  U(  and  uu  are  functions  of  the  sensor  outputs 
from  the  plant.  Like  in  Ref.  [5],  we  consider  three  sensor 
planes,  located  at  z=zsic,  k=  1,2,3.  The  sensor  planes  mea¬ 
sure  temperatures 

zsk(x,y,t)  =  0(x,y,zsbt),  k=  1,2,3.  (11) 

The  spectral  decomposition  of  the  three-dimensional  ve¬ 
locity  and  temperature  field  variables  and  the  fractional  step, 
time- splitting  method  for  integrating  the  fully  nonlinear  OB 
system  of  equations  are  described  in  detail  in  Ref.  [6].  De¬ 
tails  are  provided  in  the  references  therein,  in  particular,  by 
Marcus  [7].  Here,  only  the  addition  of  the  conductive  walls 
into  the  existing  numerical  scheme,  described  by  Eqs.  (4) 
and  (5),  will  be  mentioned  with  no  repetitions. 

Unfortunately,  a  direct  incorporation  of  Eqs.  (4)  and  (5) 
using  the  existing  fractional-step  approach  [6]  is  numerically 
unstable.  So,  the  three  heat  equations  have  to  be  solved  si¬ 
multaneously.  Referring  to  Ref.  [6],  we  add  to  the  physical 
dependent  variables  two  thermal  fields  (8f,  8U),  besides  the 
existing  fields  (u  ,v  ,w  ,p ,  8).  All  the  fields  have  independent 
variables  (jc ,y,z,t),  where  each  field  has  the  same  {x,y)  co¬ 
ordinate  ranges,  but  8f  ,  8 ,  and  8„  have  different  z  coordinate 
ranges.  We  use  three  column  vectors  to  represent  the  corre¬ 
sponding  Fourier-Chebyshev  coefficients  0U)  corre¬ 

sponding  to  the  physical  field  variables  ( 8(,8 ,  8U). 

The  finite  wall  temperature  distribution  has  a  simple 
closed-form  solution  if  the  diffusive  term  in  the  diffusive 
equations  dominates  and  the  time  derivative  term  is  negli¬ 
gible  [10].  In  general,  a  numerical  solution  is  computed.  At 
the  nth  time  step  (not  fractional  step),  the  wall  temperatures 
are  expressed  as 

0f'+1)  =  b1(a,r0('i))  +  b2M^n),  ar  =  [1,—  U, -  1,  ...f, 

^;+l)=b3(pT^n))+b4u^\  [1,1, 1,1, ... F,  (i2) 

where  (j=  1,4)  are  obtained  by  inverting  the  two  matrix 
equations  for  the  wall  heat  equations.  The  implicit  Euler  and 
tau  method  are  used  here.  For  the  lower  wall  heat  equation, 
its  lower  boundary  condition  is  the  lower  wall  control  and  its 
upper  boundary  condition  can  be  either  the  interfacial  tem¬ 
perature  matching  condition  or  the  interfacial  heat  flux 
matching  condition.  Similar  conditions  for  the  boundary  con¬ 
ditions  of  the  upper  wall  heat  equation  are  developed,  al¬ 
though  we  choose  the  interfacial  temperature  matching  con¬ 
ditions  as  the  boundary  conditions  for  solving  the  wall  heat 
equations.  The  two  scalar  controls,  u(  and  uu  are  either  func¬ 
tions  of  ( x,y,t )  (in  physical  space)  or  {kx,ky,t)  (in  Fourier 
space).  In  this  way,  the  wall  temperatures  are  explicitly 
solved  in  terms  of  the  two  control  parameters  and  the  fluid 
layer  temperature  distribution. 

We  can  now  replace  the  idealized  lower  and  upper  bound¬ 
ary  conditions  [i.e.,  8[x,y,0,t)  =  Uf(x,y,t)  and  8{x,y,\,t ) 
=  uu(x,y  ,f)]  previously  used  [6]  by  the  following  conditions, 
respectively: 
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I  Ft-  i  |  fast  fourier  transform 


FIG.  2.  The  LQG  control  loop  diagram. 

( «rD  -  ^(^DbOor)  0<"+1)  =  y  (/F Db2)4n), 

( /FD  -  y  (rFDb 3)/?)  <F+1)  =  ^(arDb4)F;!),  (13) 

where  D  is  the  first  derivative  matrix  defined  for  s  =  0, 1  for 
the  Chebyshev  coefficients  (0— >  8,8  corresponds  to  0—^  DO). 
If  the  interfacial  temperature  matching  conditions  are  used 
for  the  wall  temperatures,  then  the  interfacial  heat  flux 
matching  conditions  have  to  be  used  for  the  fluid,  and  vice 
versa.  It  can  be  shown  that,  as  >°°,  the  lower  wall 

boundary  condition  approaches  the  idealized  boundary  con¬ 
dition  used  in  our  previous  studies.  A  similar  condition  holds 
for  the  upper  wall  boundary  condition.  An  actuator  lag  for 
the  controls  is  incorporated  as  an  additional  plant  state. 

B.  Modally  distributed  LQG  compensator 

The  LQG  compensator  design  is  described  in  detail  in 
Ref.  [5],  We  will  not  repeat  the  details  here.  In  brief,  each 
compensator  (denoted  as  a  transfer  function  matrix  Kmn  be¬ 
low,  m=  1 ,2, . . .  ,NX,  n=  \  ,2, ...  ,Ny,  see  Fig.  2),  is  of  the 
linear  quadratic  Gaussian  and  loop  transfer  recovery 
(LQG-LTR)  type.  The  LTR  approximation  provides  almost 
the  full-state  feedback  performance  of  ±60°  phase  margin 
and  -6  dB  to  °°  in  lower  and  upper  gain  margins,  by  allow¬ 
ing  the  observer  weighting  parameter  /3— [11],  As  in  our 
previous  design  [6],  we  use  /3=103.  We  use  the  regulator 
weighting  parameter  7=0.1.  Note  that  -y— >0  corresponds  to 
unlimited  control  authority.  In  addition  to  the  compensator 
parameters  y  and  f3.  the  compensator  has  system  model  that 
involves  the  major  physical  parameters  of  the  problem;  the 
Rayleigh  number  as  Ra  and  the  wave  number  as  (mkx ,  nky) 
(the  asterisk  denotes  nominal  values).  When  the  nonlinear 
plant  yields  a  solution  described  by  a  2D  wave  number  array 
(mkx , nky)  for  m=  1 ,2, . . .  ,NX,  and  n  =  l  ,2, ...  ,Ny,  then  the 
compensator  consists  of  N=NxNy  single-wave-number  sub- 
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FIG.  3.  (Color  online)  Stability 
margin  in  Ra  vs  k  for  Pr=7.0 
(upper  panel)  and  correspond¬ 
ing  imaginary  eigenvalue  (lower 
panel). 


compensators.  A  £th  subcompensator  consists  of  a  Kalman 
filter  Fk  and  a  optimal  regulator  Rk  (k=  1 ,2, ...  ,N).  The  three 
levels  of  temperatures  in  Fig.  2  are  denoted  by  Zh  Z2,  and  Z3. 
The  output  from  the  k\h  compensator  is  either  one  or  two 
scalar  controls,  depending  on  one  or  two  actuator  planes 
used.  The  corresponding  plant  parameters  (without  asterisk, 
with  uncertainties)  can  differ  from  the  compensator  param¬ 
eters  (with  asterisk,  at  nominal  conditions).  Figure  2  presents 
a  block  diagram  showing  the  Fourier-decomposed  nonlinear 
plant.  All  inputs  and  outputs  are  Fourier-Chebyshev  coeffi¬ 
cients.  Each  Fourier  mode  disturbance  corresponding  to 
wave  numbers  ( mkx,nky )  is  controlled  by  a  LQG  compensa¬ 
tor  Kmn  designed  at  the  nominal  wave  numbers  ( mkx,nky ). 
The  greater  the  difference  between  the  compensator  and 
plant  parameters  while  preserving  closed-loop  stability,  the 
more  robust  the  compensator,  i.e,  the  greater  the  stability 
margins.  We  denote  the  compensator  corresponding  to  ideal¬ 
ized  thermal  boundary  condition  as  Cl,  and  the  one  incorpo¬ 
rating  finite  walls  as  C3.  The  reduced-order  compensator  [8] 
is  used  here.  Since  our  focus  is  on  robustness,  we  will  not 
deal  with  the  transient  conductive  state  nor  will  the  gain¬ 
scheduling  algorithm  be  used. 

III.  RESULTS 

A.  Linear  stability  using  a  single-plane  actuator 

In  the  linear  stability  study,  a  linearized  plant  dynamical 
model  is  used.  We  assume  the  material  properties  of  the  up¬ 
per  and  lower  walls  to  be  the  same,  corresponding  to  alumi¬ 
num  at  room  temperature.  The  walls  have  thickness  1/10  of 
the  fluid  depth.  For  the  working  fluid  using  water  we  have 
Pr=7.0,  Ki=K,=  400,  and  ki=k1=670.  Using  air  we  have 
Pr=0.7,  K/=KU  =  8800,  and  Ki=ku= 4.2  instead. 


1.  Major  parameter  uncertainties 

Our  nominal  model  consists  of  one  actuator,  located  at 
outer  surface  of  the  lower  wall,  at  z=-dt  and  three  sensor 
planes,  located  at  zsl=0.2,  zj2  =  0.5,  and  zj3  =  0.8.  These  are 
optimal  sensor  locations  according  to  our  previous  analysis 
[5].  The  nominal  Ra  =20  Rac,  where  Ra,  ~  1707.762  is  the 
critical  Rayleigh  number.  The  LQG  compensator  is  designed 
at  the  nominal  values  (denoted  by  an  asterisk).  In  the  non¬ 
linear  simulations  [6-8]  the  compensator  consists  of  an  array 
of  single-/:  controllers.  The  nominal  wave  number  k  is  an 
array  evenly  spaced  wave  numbers  covering  the  entire  un¬ 
stable  band  at  Ra*.  In  the  linear  stability  analysis,  the  nomi¬ 
nal  wave  number  k*  is  a  prescribed  parameter  that  can  be 
varied  across  the  unstable  band.  Here  the  compensator  is 
assumed  to  be  designed  at  (k  ,  Ra).  At  nominal  condition, 
the  plant  model  has  k=k*  and  Ra=  Ra  .  For  robustness  inves¬ 
tigation,  we  vary  k  and  Ra  from  their  nominal  values. 

The  eigenvalue  analysis  indicates  that  the  closed-loop 
system  using  the  steady-state  LQG  compensator  is  always 
stable  at  the  nominal  condition,  i.e.,  at  k=l<  and  Ra=Ra 
even  when  Ra  is  as  large  as  100Rac.  However,  at  high  Ra 
exceeding,  say,  20Rar,  the  least  stable  eigenvalue  becomes 
sensitive  to  small  changes  in  the  prescribed  parameter.  This 
is  an  indication  of  weakening  of  robustness  at  large  Ra. 

In  the  following,  we  vary  the  plant  parameter  Ra  from  the 
nominal  value  Ra  =20Raf  to  see  how  far  the  system  can 
tolerate  the  change  before  becoming  unstable.  The  result  pre¬ 
sents  the  stability  margin  ARa=  Ra-Ra  in  terms  of  the  wave 
number  k* .  In  doing  so  we  assume  that  the  plant  wave  num¬ 
ber  k  is  unchanged  from  the  nominal  value  k*.  In  Figs.  3  and 
4,  we  show  the  results  Ra  vs  k  for  the  case  at  Pr=7.0  and 
Pr=0.7,  respectively.  The  margin  ARa  is  the  range  above  20. 
For  the  two  values  of  Pr,  the  results  appear  very  similar.  The 
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Wave  number 


FIG.  4.  (Color  online)  Stability 
margin  in  Ra  vs  k  for  Pr=0.7  (up¬ 
per  panel)  and  corresponding 
imaginary  eigenvalue  (lower 
panel). 


margin  in  Ra  (upper  panels)  dips  to  a  minimum  near  /c=  13. 
Although  the  minimal  margin  occurs  at  approximately  four 
times  the  critical  wave  number  kc~3.117,  in  the  nonlinear 
simulations  presented  later,  we  see  no  sign  that  short-wave 
modes  are  being  excited.  In  fact,  nonlinearity  is  stabilizing 
by  cascading  disturbances  to  shorter  and  shorter  harmonic 
modes  that  promote  diffusions.  Beyond  k=  13,  the  modal  dif¬ 
fusive  terms  become  very  large  and  stabilizing;  the  Ra  mar¬ 
gin  is  expected  to  increase  unboundedly.  The  pair  of  imagi¬ 
nary  eigenvalues  with  zero  real  part  is  shown  in  the  lower 
panels.  The  result  shows  that,  for  positive  plant  uncertainty 
in  Ra  (i.e.,  Ra>  Ra  ),  the  closed-loop  behavior  first  becomes 
unstable  to  an  oscillatory  mode  at  lower  wave  number  and  to 
a  monotonic  mode  at  higher  wave  number.  For  negative 
plant  uncertainty  in  Ra  (i.e.,  Ra<  Ra  ),  the  closed-loop  be¬ 
havior  is  always  stable. 

We  turn  to  the  wave  number  margins.  Here  we  assume 
there  is  no  uncertainty  in  Ra,  i.e.,  Ra=Ra  =20Rac.  We  vary 
the  nominal  wave  number  k*  and  at  each  value  of  k  we 
determine  the  Ak  margins,  A k=k-k*,  where  there  are  two 
values  of  Ak,  the  upper  and  lower  Ak  stability  limit.  Again, 
we  consider  the  cases  of  Pr=7.0  (Fig.  5)  and  Pr=0.7  (Fig.  6). 
Again,  we  see  that  the  behavior  for  both  cases  are  very  simi¬ 
lar.  In  both  cases,  the  upper  Ak  margin  (upper  curve)  in¬ 
creases  steeply  beyond  k  =  8.  The  lower  k  margin  (lower 
curve)  behaves  in  a  rather  complicated  fashion.  The  margin 
diverges  between  k* =  4.0  and  5.0  but  between  k* = 5.0  and 
the  Ak  minimum  near  k  =  13,  the  margin  is  quite  flat.  Then, 
beyond  k ’  =  13,  the  lower  Ak  margin  increases  linearly.  It 
appears  that  the  stability  margin  is  weakest  near  k  =  13, 
which  means  that  the  closed-loop  system  is  least  robust  with 
respect  to,  or  most  vulnerable  to,  unstable  mode  onset  for  a 
spatial  length  scale  corresponding  to  k  =  13,  or  a  wavelength 
X  =2tt Ik" ,  approximately  1/2  of  the  layer  depth. 


We  emphasize  that  the  three-level  sensor  configuration  is 
important  for  robustness.  Decreasing  the  number  of  sensor 
planes  or  moving  the  outer  sensor  planes  towards  the  wall 
will  erode  the  state  observability.  As  a  result,  closed-loop 
stability  is  only  achievable  at  much  lower  Ra*. 

2.  Actuator  lag 

A  potential  destabilizing  effect  arises  from  the  actuator 
lag.  Both  finite  conductivity  and  electronic  processing  in  the 
thermal  actuator  contribute  to  the  actuator  lag.  A  first-order 
lag  effect  is  incorporated  in  the  linearized  plant  model  and 
later  in  the  nonlinear  plant,  described  by  the  transfer  function 
(TF)  Gd=utd/(s  +  cod)  where  s  is  the  Laplace-transformed 
value.  This  TF  has  to  be  modeled  in  between  the  compensa¬ 
tor  output  ua  and  the  plant  input  Because  of  the  lag,  we 
have 


ui  +  wdui  =  o)du0.  (14) 

For  large  wd,  or  small  lag  time  constant  tc1=2ttI  a>d,  the  effect 
is  small,  i.e.,  But  for  large  rd,  the  lag  can  be  desta¬ 

bilizing.  To  model  the  actuator  lag,  we  write  the  linear  plant 
state-space  equation  as  x=Ax+Bm,-,  where  u,  is  given  in  Eq. 
(14).  The  compensator  output  ua  can  be  used  as  the  input  to 
the  plant  model,  modified  as 


X 

= 

"a 

_o 

B 

X 

+ 

0 

Uj 

Uj 

y  =  [c 


(15) 
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FIG.  5.  (Color  online)  Margins 
A k  vs  k*  for  Ra=20Rac,  Pr=7.0. 


In  Fig.  7,  we  show  the  destabilizing  effect  of  the  actuator  lag 
by  plotting  the  real  part  of  the  least  stable  eigenvalue 
(growth  rate)  versus  wave  number.  In  this  case,  both  plant 
and  compensator  parameters  are  set  to  the  same  condition, 
with  k=k* ,  Ra=Ra  =20Raf.  The  left  panel  shows  the  case 
for  Pr=7.0  and  the  right  panel  shows  the  case  for  Pr=0.7.  In 


both  cases,  the  least-stable  eigenvalue  (real)  is  shown  as  a 
function  of  the  wave  number  (A:)  for  three  lag  time  constants, 
7^=0.001,  0.003,  and  0.005.  These  are  in  nondimensional 
time  unit.  In  the  wave-number  band  between  k=3  and  k=  8 
for  the  case  of  Pr=7.0,  a  sufficiently  large  lag  time-constant 
causes  instability.  For  the  case  Pr=0.7,  however,  the  lag  has 


FIG.  6.  (Color  online)  Margins 
A k  vs  k*  for  Ra=20Rar,  Pr=0.7. 
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(a)  Pr=7.0 


(b)  Pr=0.7 


FIG.  7.  (Color  online)  Least 
stable  eigenvalue  vs  k  with  actua¬ 
tor  lag. 


no  destabilizing  effect.  For  the  case  Pr=0.7,  rd=  0.003  al¬ 
most  drives  the  Fourier  mode  at  k  =  5  unstable.  For  the  case 
Pr=0.7,  the  actuator  lag  is  more  destabilizing  for  higher 
wave  number.  The  eigenvalue  peaks  near  £=13. 

The  nondimensional  time  scale  is  d2/ k.  For  water  and 
layer  thickness  <7=0.5  cm,  one  nondimensional  time  unit 
=  176  s  at  standard  temperature  pressure  (STP).  For  air  k  is 
about  145  times  that  of  water.  In  order  to  get  the  same  non- 
dimensional  time  scale  d  for  air  has  to  be  12  times  that  of  d 
for  water.  For  the  nondimensional  time  scale,  if  the  layer 
thickness  is  0.5  cm,  say  for  water  at  STP,  a  physical  lag  rd 
=  0.001  corresponds  to  0.18  s,  roughly;  and  that  for  air  is 
only  0.001  s.  One  can  always  stretch  the  physical  lag  time 
for  the  same  value  of  rd,  by  increasing  the  layer  depth  of  air. 
In  order  to  maintain  the  same  Ra  for  air  (Ra  is  proportional 
to  d2/  k  for  a  given  temperature  difference),  however,  this 
implies  decreasing  the  temperature  difference  for  air. 


3.  Sensor  plane  depth  uncertainty 

The  sensor-plane  depth  uncertainty  turns  out  to  be  a  dif¬ 
ficult  problem.  The  nominal  sensor  plane  depths  are  z\=0.2, 
z2  =  0.5,  and  z3  =  0.8  (scaled  by  the  fluid  layer  thickness). 
Consider  the  compensator  model  to  have  these  nominal  val¬ 
ues  but  the  plant  model  has  corresponding  plane  depths  of 
Zi=z\+tSz\,  z2  =  z2  +  Az2,  and  z3=73  +  Az3  where  the  perturba¬ 
tions  are  due  to  uncertainties.  We  assume  each  of  Az,  is  a 
zero-mean  Gaussian  number  with  standard  deviation  0.01. 
Consider  the  case  Pr=7.0,  Ra=Ra  =2000,  £=£*  =  3.5.  We 
use  the  C3  compensator  for  the  three-layer  plant  model.  Fig¬ 
ure  8  shows  the  results  for  5000  Monte  Carlo  runs.  For  each 
realization  of  error  (Azi , Az2, Az3),  the  closed-loop,  least 
stable  eigenvalue  is  computed.  If  the  real  part  is  less  than  or 
equal  to  zero,  a  dot  is  shown  in  the  plot.  The  ensemble  shows 
where  the  stable  region  lies  on  the  3D  error  space.  Panels 


FIG.  8.  (Color  online)  Stable 
region  in  the  3D  error  space  from 
5000  Monte  Carlo  runs. 
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Wave  number 


FIG.  9.  (Color  online)  Stability 
margin  in  Ra  vs  k  for  Pr=7.0  (up¬ 
per  panel),  and  corresponding 
imaginary  eigenvalue  (lower 
panel). 


(a)-(c)  are  projection  of  the  points  (stable)  onto  the  Azt  vs 
A z2  plane,  Az2  vs  Az3  plane,  and  A Zi  vs  Az;,  plane,  respec¬ 
tively.  The  last  panel  shows  the  3D  plot  of  the  points.  The 
unstable  points  are  not  shown.  The  figures  use  a  uniform 
scale  per  axis  from  -0.04  to  0.04.  The  results  indicate  that  (i) 
the  stable  region  in  the  3D  error  space  is  quite  small;  (ii)  A z3 
appears  to  be  the  least  significant  error  among  the  three, 
relatively  speaking.  It  means  the  sensor  plane  farthest  from 
the  actuator  plane  can  tolerate  largest  uncertainty;  (iii)  both 
Azi  and  Az2  are  important.  However,  the  shape  of  the  stable 
region  is  not  symmetric  with  respect  to  the  sign  of  error.  It 
appears  less  destabilizing  if  the  first  and  second  sensor 
planes  are  moving  away  from  the  actuator,  more  destabiliz¬ 
ing  if  both  are  moving  towards  the  actuator.  It  is  worse  if  one 
is  moving  away  but  the  other  is  moving  towards  the  actuator. 
The  sensor-plane  depth  errors  can  be  a  significant  challenge 
in  the  laboratory  implementation. 

B.  Linear  stability  using  two  actuator  planes 

It  is  of  interest  to  see  if  there  will  be  improvement  in  the 
robustness  by  the  addition  of  an  extra  actuator  plane  on  the 
outer  face  of  the  upper  wall.  This  addition  makes  the  plant  a 
three-input-two-output  system.  We  choose  to  show  the  case 
Pr=7.0  only.  In  Fig.  9,  we  show  the  Ra  margin  in  k  as  in  the 
nominal  case.  The  margin  curves  are  very  similar  to  the  ones 
before.  Only  a  very  slight  improvement  in  margins  is  evi¬ 
dent.  The  improvement  is  not  significant.  Compared  the  sec¬ 
ond  panel  between  Figs.  3  and  9,  the  k  band  for  imaginary 
part  of  the  least-stable  eigenvalue  becomes  significantly 
smaller.  In  Fig.  10,  we  show  the  upper  and  lower  k  margins 
vs  k.  Again,  the  improvement  with  an  additional  actuator  is 
small. 


1.  Singular-value  bounds 

Doyle  and  Stein  [11]  develop  an  approach  to  use  the  sin¬ 
gular  values  bounds  for  relative  stability  measure  in  finite¬ 
dimensional,  linear-time-invariant  (FDLT)  systems.  The  con¬ 
ditions  give  bounds  that  guarantee  stability  but  these  are  not 
necessarily  tight  bounds,  therefore,  more  conservative.  To 
illustrate  the  idea,  it  is  more  convenient  to  use  the  transfer 
function  (TF)  notation.  The  TF  between  input  and  output  of 
the  state-space  system  (A,B,C,D),  corresponding  to  the 
standard  dynamical  and  output  equations  x=Ax+Bu  and  y 
=  Cx+Du  (after  Laplace  transform)  is  v(.v)  =  G(,s]u(.s],  where 
G(s)  =  C(sI-A)-‘B+D.  Here,  we  use  G(,s]  to  denote  the 
plant  TF.  In  the  following,  both  plant  and  compensator  TFs 
correspond  to  a  single-wave-number  model.  There  are  two 
common  ways  to  model  uncertainties.  One  is  by  additive 
perturbation  to  the  nominal  G(.v),  G'=G  +  AG.  The  other  is 
by  multiplicative  perturbation  to  the  nominal  G(.v),  G'  =  (I 
+L)G,  where  L  is  a  multiplicative  error  model  derived  based 
on  G'  and  G.  Here,  the  latter  way  is  adopted.  To  keep  the 
uncertainties  within  bounds,  we  require  7r(L)  <  f„,(oj)  for 
some  prescribed  function  of  frequency,  Here  a  and  cr 

denote  the  upper  and  lower  bound  singular  values,  respec¬ 
tively.  The  LQG  compensator,  denoted  by  K(.v)  is  defined  by 
state  space  system  ((A*-KyC*-B*Kc) ,Ky,-KeD*),  where 
K f  and  K  are  the  filter  and  regulator  gains,  respectively. 
Otherwise,  the  matrices  are  identical  to  those  of  the  linear¬ 
ized  plant. 

There  are  two  ways  consider  breaking  the  loop  of  a  com¬ 
pensated  plant.  One  is  by  breaking  the  loop  at  the  plant  out¬ 
put  [see  Fig.  11(a)],  This  case  gives  a  compensated  plant  TF 
equal  to  GK(s)  and  the  other  is  by  breaking  the  loop  at  the 
plant  input  [see  Fig.  11(b)].  In  both  figures,  the  dashed  line  in 
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FIG.  10.  (Color  online)  Mar¬ 
gins  A k  vs  k*  for  Ra=20Rac,  Pr 
=  7.0. 


the  loop  denote  the  broken  loop.  This  case  gives  a  compen¬ 
sator  plant  TF  of  KG(s)  (note  that  the  matrix  columns  cor¬ 
respond  to  the  inputs  and  matrix  rows  correspond  to  the  out¬ 
puts).  In  Ref.  [5],  we  considered  a  plant  with  one  actuator 
plane  (one  input)  and  three  sensor  planes  (three  outputs). 
Therefore  the  dimensions  of  G(.v)  is  3  X  1 .  The  compensator 
takes  the  three  plant  outputs  as  measurements  to  produce  one 
control  therefore  the  dimension  of  K(s)  is  1X3.  The  advan¬ 
tage  by  breaking  the  loop  at  the  plant  input  is  that  we  obtain 
a  single-input-single-output  (SISO)  GK(s).  Therefore,  as 
demonstrated  in  Ref.  [5],  the  relative  stability  can  be  effec¬ 
tively  analyzed  using  gain  and  phase  margins,  based  on  the 
classical  Nyquist  criterion. 

Here,  the  Doyle  and  Stein  condition  provides  an  indepen¬ 
dent  means  to  assess  the  stability  margins,  apart  from  the 
parameter  margin  curves  from  the  direct  closed-loop  compu¬ 
tation.  We  consider  two  actuator  planes  and  three  sensor 
planes.  The  Nyquist  criterion  can  no  longer  be  applied.  Con¬ 
sider  breaking  the  loop  at  the  plant  output  (the  argument  is 


G 

G 

u 

|  z  u| 

K 

K 

(a)  Breaking  the  Loop  at  (b)  Breaking  the  Loop  at 

Plant  Output  tf=gk  Plant  input ,  tf=  kg 

FIG.  1 1 .  Block  diagram  showing  breaking  the  compensated 
plant  loop. 


equally  valid  for  breaking  the  loop  at  the  plant  input.  In  the 
present  case,  however,  the  plant  has  two  plant  inputs  versus 
three  outputs,  we  caution  that  the  3X3  matrix  GK  is  rank 
deficient  by  1).  Doyle  and  Stein  [11]  noticed  that  to  guaran¬ 
tee  stability  in  the  presence  of  uncertainties,  G'K(s)  has  to 
preserve  the  system  dimensions  of  GK(s).  If  the  compen¬ 
sated  plant  has  equal  inputs  and  outputs,  then  the  generalized 
Nyquist  criterion  on  the  determinant  is  that  the  number  of 
encirclement  of  det(I+G'K)  remains  unchanged  from 
det(I  +  GK).  Preserving  the  system  dimensions  means  that 
the  lower  bound  of  the  singular  value  of  I  +  G '  K  is  positive 
definite.  That  is, 

0  <  tr[I  +  (I  +  L)GK].  (16) 

Since  (I+GK)  is  full  rank,  we  pull  it  out  as  the  common 
factor.  Then,  we  require 

0  <  of  I  +  LGK(I  +  GK)’1]  (17) 

for  all  to  and  L,  or 

ct[GK(I  +  GK)-1]  <  l/€m.  (18) 

Using  the  matrix  identity 

GK(I  +  GK)-1  =  [I  +  (GK)-1]-1  (19) 

the  stability  condition  according  to  Ref.  [11]  is 

otL(®)]  <  pfl  +  (GK)’1],  (20) 

where  the  right-hand  term  is  referred  to  as  the  matrix  return 
difference  of  the  loop  and  €m(<u)  represents  the  singular- 
value  bound  of  the  error  model  <t(L(oj)).  Note  that  GK  com¬ 
mutes  with  (I  +  GK)-1.  In  the  following,  we  apply  the  Doyle- 
Stein’s  stability  condition  to  investigate  two  separate 
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FIG.  12.  (Color  online)  Singular-value 
bounds  with  the  actuator  error  model. 


problems.  One  is  the  actuator  lag.  The  other  is  the  finite  wall 
thermal  properties.  The  actuator  is  based  on  a  second-order 
lag  model  according  to  Ref.  [12].  Let  L  (for  this  case  a  sca¬ 
lar)  be  given  by 


a-(Uaij)  = 


s 2  +  2  £;coas 
s2  +  2  £was  +  co2 


S=j(0. 


(21) 


A  damping  factor  £=0.15  is  used  but  £  is  not  a  sensitive 
parameter.  Decreasing  £  to  0.05  gives  roughly  the  same  re¬ 
sults.  In  Figs.  12(a)  and  12(b),  we  show  the  case  <wn  =  5000 
and  in  Figs.  12(c)  and  12(d),  we  show  oj„  =  2000.  These  are 
nondimensional  frequencies.  The  corresponding  physical  fre¬ 
quencies  are  typically  two  orders  of  magnitude  lower.  The 
left  panel  corresponds  to  k  =  6.5  and  the  right  panel  to  k 
=  13.0.  The  higher  wave-number  system  appears  more  robust 
with  respect  to  the  actuator  lag  than  the  lower  wave  number. 
In  each  panel  the  two  dashed  curves  correspond  to  the  lower- 
bound  singular  values  [see  the  right-hand  side  of  inequality 
(20)].  The  solid  curve  shows  <f(L).  The  stability  margin  cor¬ 
responds  to  the  gap  between  the  lower  dashed  curve  and  the 
solid  curve.  The  larger  the  gap,  the  greater  the  margin.  The 
result  appears  fairly  consistent  with  those  of  Fig.  7.  The 
upper-bound  singular  value  of  the  actuator  lag  is  shown  in 
the  lowest  curve  in  each  panel  of  the  figure.  Now,  the  lag 
time  constant  is  tci  =  2it/  u>a.  The  smaller  oj„  the  greater  the 
lag  and  the  more  destabilizing  the  lag  effect  becomes.  The 
closed-loop  system  remains  stable.  The  physical  implication 
is  this.  For  water  as  fluid,  an  actuator  lag  time  constant  of 
0.5  s  will  not  trigger  instability.  For  air  at  STP,  this  time 
constant  is  shortened  by  a  factor  of  roughly  150  times.  It 
appears  that  if  the  physical  dimensions  are  unchanged,  con¬ 
trolling  convection  in  air  requires  a  much  higher-bandwidth 
controller  than  in  water.  Next,  we  address  the  problem  of 
uncertainties  of  the  finite  wall  thermal  properties.  Here,  we 


construct  a  plant  error  model,  denoted  by  AA  and  AB  (where 
C  and  D  are  not  affected)  by  the  following.  Let  Aku=Ak( 
=  A k  and  AKII  =  AK(  =  AK.  We  compute  the  derivatives 
dA/dK,  i?B  /  dK,  dA  I  dK,  and  dH/dK.  So 

dA  dA 

AA  =  — A  k+  — A  K, 
dK  dK 


AB  =  — A  k+ — A  K.  (22) 

dK  dK 

Note  that  the  error  model  affects  only  the  plant  and  the  com¬ 
pensator  remains  at  nominal  parameters.  From  the  error 
model,  it  is  somewhat  tedious,  but  quite  straight  forward  to 
generate  the  TF  version  of  the  error  L  by  keeping  the  first- 
order  error  terms  (the  derivation  will  not  be  produced  here). 
Figures  13(a)  and  13(b)  shows  the  singular-value  bounds  for 
the  same  wave  numbers  k=  6.5  and  13.0,  respectively.  In  this 
case,  we  let  A k  and  A K  to  be  25%  of  the  nominal  values. 
The  dashed  curves  are  the  same  as  in  Fig.  12.  The  pair  of 
lower  solid  curves  represent  the  two  singular  values  of  L(oj). 
The  solid  and  dashed  curves  very  barely  intersect  for  the  k 
=  6.5  case  near  oj~  700.  It  appears  that  stability  is  preserved 
in  the  presence  of  25%  uncertainties  in  the  wall  thermal 
properties  (both  diffusivity  and  conductivity),  about  the 
nominal  values. 


C.  3D  nonlinear  closed-loop  simulations 
1.  Model  setup  and  assumptions 

The  compensator  is  linear,  set  at  the  nominal  parameters 
and  the  plant  model  is  3D,  nonlinear  with  uncertainties  about 
the  nominal  parameters.  The  closed-loop  system  in  the  pres¬ 
ence  of  uncertainties  is  investigated  by  examining  the 
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(a) 


(b) 


FIG.  13.  (Color  online)  Singular- value  bounds  with  the  finite 
wall  property  error  model. 

closed-loop  time  response  of  the  upper  and  lower-wall 
Nusselt  numbers.  Like  in  Refs.  [6-8],  the  measure  of  the 
residual  convection  can  be  effectively  shown  using  the  Nus¬ 
selt  number  plots.  The  Nusselt  numbers  at  the  upper  and 
lower  wall  measure  the  ratios  of  total  heat  flux  (convective 
and  conductive)  to  the  conductive  heat  flux  leaving  and  en¬ 
tering  the  fluid  layer.  As  Nusselt  number  approaches  the 
value  1.0,  the  fluid  layer  convection  is  removed.  Robustness 
is  determined  based  on  how  much  mismatches  can  be  toler¬ 
ated  before  the  closed-loop  system  is  unstable. 

Although  the  parameter  field  is  extremely  vast,  we  have 
simulated  a  large  number  of  cases  but  for  limited  space  only 
a  selective  number  of  representative  cases  are  presented. 
Since  the  focus  of  this  paper  is  the  robustness  of  the  closed- 
loop  system  with  respect  to  parameter  mismatches,  the 
Rayleigh  number  is  set  at  a  constant  during  the  simulation. 
No  gain-schedule  algorithm  [8]  will  be  engaged.  The  nomi¬ 
nal  parameters  are  set  at  Pr=7.0,  Ra  =2X104,  k*x  =  k\~  1.0. 
Unlike  the  linear  analysis,  these  are  fundamental  wave  num¬ 
bers.  The  nominal  geometric  and  material  properties,  sensor 
and  actuator  configurations  are  the  same  as  in  the  linear 
analysis  (see  first  paragraph  of  Sec.  Ill  A).  The  simulations 
correspond  to  32X32  Fourier  modes  (horizontal).  The  inte¬ 
gration  and  output  sampled  step  is  A/ =  0.004.  The  simulation 
period  T  is  0.4  time  units.  Only  one  set  of  initial  condition  is 
used,  corresponding  to  the  residual  state  of  the  closed-loop 
simulation  at  1=0.4  presented  in  Ref.  [8].  We  use  two  ver¬ 


sions  of  the  LQG  compensator.  Cl  and  C3,  designs  based  on 
a  one-layer  and  three-layer  model,  respectively.  Both  are 
based  on  a  reduced-order  linear  model  consisting  of  only 
eight  vertical  complex  modes  (contrast  to  64  vertical  real 
modes  in  the  plant  model). 

2.  Finite  wall  properties 

We  selected  a  handful  of  cases  to  characterize  the  closed- 
loop  behavior  (both  stable  and  unstable  cases  included).  Fig¬ 
ure  14(a)  shows  the  upper  and  lower  Nusselt  number  re¬ 
sponses  when  the  compensator  Cl  is  used.  The  plant 
parameters  are  set  at  nominal.  Note  that  the  compensator  has 
idealized  thermal  boundary  conditions  whereas  the  plant  has 
finite  walls  incorporated.  The  upper  (dashed)  and  lower 
Nusselt  number  (solid)  show  convection  damped  out  in  time. 
We  then  switch  to  compensator  C3,  the  time  response  (not 
shown  here)  is  visibly  indistinguishable  from  that  of  Fig. 
14(a). 

Keeping  the  plant  at  nominal  condition,  now  we  reduce 
the  upper  and  lower  wall’s  thermal  conductivity  and  diffu- 
sivity  values  each  by  50%.  The  closed-loop  system  is  stable 
with  Cl  as  well  as  C3. 

Next,  we  increase  both  wall  thicknesses  from  the  nominal 
value  to  df=du  =  0.15.  Compensator  Cl  is  not  capable  of  sup¬ 
pressing  convection  any  longer.  Figure  14(b)  shows  the  ini¬ 
tial  response.  Eventually,  the  solution  diverges.  Using  com¬ 
pensator  C3,  with  the  nominal  wall  thicknesses  at  0.1,  C3 
still  cannot  damp  out  convection.  In  this  case,  the  solution 
takes  longer  to  diverge  (plot  not  shown). 

We  conclude  that  the  mismatches  in  wall  thermal  proper¬ 
ties  between  the  compensator  and  plant  up  to  50%  is  easily 
tolerated.  But  the  mismatch  in  wall  thickness  is  significantly 
more  sensitive. 


3.  Simulations  with  actuator  lag 

In  evaluating  the  actuator  lag,  we  reset  all  the  plant  pa¬ 
rameters  to  nominal  values.  To  incorporate  the  first-order  ac¬ 
tuator  lag,  the  last  of  Eqs.  (15)  has  to  be  incorporated  into  the 
time-splitting  algorithm  of  the  3D  nonlinear  plant  model.  We 
rewrite  the  equation  for  the  lower  plane  actuator  (similar 
equation  for  the  upper  plane  actuator), 


-  -  Mdufi  +  Mduto- 


(23) 


Subscript  o  denotes  the  output  from  the  compensator  and  i 
denotes  the  input  to  the  plant;  the  lag  frequency  tod=2 Trlrd. 
We  attempted  to  integrate  the  above  first-order  equation  by 
the  explicit  Euler  scheme.  The  closed-loop  solution  is  nu¬ 
merically  unstable.  We  switch  to  the  implicit  Euler  scheme. 


(n+l)  _  llt "i  +  ^t(°dueo 


(n+1) 


1  +  A  t(ori 


(24) 


which  is  numerically  stable.  For  a  lag  time  constant  r„  as 
large  as  3  X  10-2,  Cl  successfully  stabilizes  the  system.  Sur¬ 
prisingly,  this  time  constant  is  significantly  higher  than  that 
predicted  by  the  linear  stability  result  (see  Fig.  7).  With  C3, 
the  performance  is  only  slightly  better.  We  increase  the  lag 
time  constant  to  rn=0.1,  ten  times  as  large.  In  Fig.  14(c),  we 
show  the  time  response.  The  closed-loop  system  is  unstable. 
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(a) 


FIG.  14.  (Color  online)  Case 
studies  in  the  nonlinear 
simulations. 


4.  Major  parameter  uncertainties 

The  most  benign  major  parameter  mismatch  by  far  is  the 
Prandtl  number.  At  nominal  conditions,  we  can  reduce  Pr 
from  7.0  to  0.7  in  the  plant  (compensator  remains  at  Pr 
=  7.0)  without  destabilizing  the  closed-loop  time  response. 

For  Rayleigh  number,  previous  results  indicate  that  the 
higher  the  nominal  value  Ra*  (here  Ra  =2X104),  the 


smaller  the  plant  uncertainty  in  Ra  can  be  tolerated.  Here,  we 
let  Ra=  1 ,05Ra  (5%  uncertainty).  Cl  is  capable  of  stabiliz¬ 
ing  the  no-motion  state.  However,  when  we  set  Ra=  I .  I  Rae, 
now  Cl  is  fighting  very  hard.  We  show  the  time  response  in 
Fig.  15(a).  Considerable  improvement  is  demonstrated  by 
using  C3  instead,  as  the  time  response  shown  in  panel  (b). 
But  the  system  is  still  unstable.  The  Ra  margin  in  the  3D 


FIG.  15.  (Color  online)  More 
case  studies  in  the  nonlinear 
simulations. 
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nonlinear  model  appears  weaker  than  that  predicted  by  the 
linear  stability  model  (see  Fig.  3). 

5.  Simulations  with  sensor  plane  level  mismatch 

Last,  we  consider  the  sensor  plane  level  mismatch.  The 
sensor  level  appears  to  be  a  sensitive  parameter  to  the  closed- 
loop  stability.  Even  though  we  cannot  run  as  many  Monte 
Carlo  cases  as  in  the  linear  stability  study  (see  Fig.  8),  we  ran 
several  cases  with  level  mismatch  with  standard  deviation 
(STD)  ±0.02  randomly  (zero  mean,  Gaussian)  added  to  the 
three  nominal  levels.  The  closed-loop  responses  remain 
stable.  Flowever,  when  we  increase  the  mismatch  STD  to 
±0.05,  the  closed-loop  system  is  unstable.  Figure  15(c) 
shows  the  time  response  for  a  case,  with  zs  =  0.15,0.55,  and 
0.85). 

In  examining  the  time  responses  of  the  Nusselt  numbers, 
it  seems  puzzling  at  first  to  see  that  the  convective  distur¬ 
bance  is  damped  out  to  a  very  small  amplitude  in  all  cases 
initially,  but  the  control  action  cannot  sustain  stability  in 
some  cases.  The  explanation  is  as  follows.  The  initial  time 
responses  depend  on  the  initial  condition  of  the  states.  Since 
only  one  set  of  initial  states  is  used,  it  is  not  surprising  that 
the  initial  responses  for  all  cases  are  similar.  The  simulation 
time  appears  adequate  for  the  unstable  modes  to  re-organize. 
The  main  point  is  that  we  have  to  simulate  long  enough  to 
pass  the  initial  transient  period.  The  asymptotic  response  is 
what  determine  stability. 

There  is  no  specific  mention  about  how  the  sensors  can  be 
implemented.  It  should  be  understood  that  the  conventional, 
invasive-type  of  temperature  sensors  are  probably  not  prac¬ 
tical.  Through  private  communications,  however,  it  comes  to 
our  knowledge  that  certain  infrared  (IR)  optical  temperature 
sensing  techniques  (remote  sensing)  are  available.  Such  op¬ 
tical  remote-sensing  method  can  probe  temperatures  at  vari¬ 
ous  depths  in  the  layer,  to  very  high  accuracy.  If  the  field-of- 
view  of  the  sensor  cannot  cover  the  total  horizontal  span  of 
the  entire  layer,  a  high-frequency  scanning  technique  can  be 
used  to  reconstruct  the  temperature  field.  For  laboratory 
implementations,  the  remote- sensing  method  in  measuring 
the  temperatures  should  be  further  investigated. 

IV.  CONCLUSION 

This  study  is  to  assess  the  amount  of  mismatches  in  the 
plant  parameters  (for  a  particular  sensor  and  actuator  con¬ 
figuration)  that  can  be  tolerated  by  the  LQG  compensator, 
before  the  closed-loop  system  turns  unstable.  This  assess¬ 
ment  is  an  important  step  towards  any  potential  future  labo¬ 
ratory  implementation.  The  assessment  is  done  by  keeping 
the  compensator  operating  at  the  nominal  values  and  intro¬ 
ducing  mismatches  to  the  plant  model. 

Based  on  the  results  from  both  the  linear  stability  study 
and  the  nonlinear  time-domain  simulations,  we  draw  the  fol¬ 
lowing  conclusions:  (i)  Introducing  the  finite  wall  to  the  fluid 
layer  does  not  have  a  significant  impact  in  altering  the 
closed-loop  stability  properties.  The  thermal  conductive  and 
diffusive  properties  of  the  walls  in  the  study  correspond  to  a 
good  conducting  material.  The  case  of  poor  conducting  ma¬ 
terial  has  not  been  considered.  Therefore  the  idealized  ther¬ 


PHYSICAL  REVIEW  E  73,  046307  (2006) 

mal  boundary  condition  used  in  Refs.  [5-8]  appears  ad¬ 
equate.  The  current  results  indicate  that  using  the 
compensator  Cl  (with  idealized  boundary  conditions)  versus 
the  finite-wall  compensator  C3  does  not  make  any  significant 
difference  at  all.  For  potential  laboratory  implementation,  the 
dynamical  model  with  idealized  boundary  conditions  is  prob¬ 
ably  adequate,  (ii)  Using  two  actuator  planes  on  both  walls 
does  not  show  significant  improvement  in  performance  over 
using  one  actuator  plane  on  the  lower  wall.  Using  two  actua¬ 
tor  planes  reduces  the  load  carried  by  one  actuator  plane,  (iii) 
For  the  parameter  mismatches  in  Rayleigh  number,  wall  con¬ 
ductive  and  diffusive  properties,  wall  thickness  and  sensor 
level  locations,  the  linear  stability  results  shows  considerably 
larger  margins  than  those  obtained  from  the  fully  nonlinear 
simulations.  For  the  actuator  lag,  however,  it  is  the  other 
way.  The  nonlinear  results  show  more  margins  that  from  the 
linear  results.  Both  linear  and  nonlinear  results  indicate  that 
the  Pr  mismatch  is  insignificant.  Mismatched  Pr  value  any¬ 
where  in  the  range  of  0. 7-7.0  is  tolerated,  (iv)  The  more 
sensitive  mismatches  arise  from  the  wall  thicknesses  and  the 
sensor  level  locations.  Both  linear  analysis  and  nonlinear 
simulations  indicate  that  only  small  mismatches  in  these  pa¬ 
rameters  can  be  tolerated  at  Ra=11.7Rac  (nominal  value). 
Reducing  the  nominal  Ra  will  reduce  the  sensitivity  for 
given  thickness  and  level  mismatches. 

In  conclusion,  the  LQG  compensator  design  is  adequate 
to  suppress  convection  in  the  vicinity  of  10-15  times  the 
critical  Ra.  The  implementation  of  the  actuator  does  not 
seem  to  be  a  challenge.  The  study  shows  that  significant 
actuator  lag  can  be  tolerated.  The  major  challenge  appears  to 
be  in  the  sensors.  The  problem  of  sensitivity  of  stability  mar¬ 
gins  to  the  sensor  levels  may  require  a  novel  temperature 
measurement  technique  to  resolve.  Rather  than  the  invasive 
method,  it  is  desirable  to  use  the  IR  remote- sensing  tech¬ 
nique,  coupled  with  a  scanner  approach.  In  principle,  such 
technique  can  deliver  temperature  measurements  at  several 
levels  of  the  fluid  layer  simultaneously,  rapidly  and  accu¬ 
rately.  Further  investigation  of  the  advanced  sensing  method 
is  necessary  for  potential  laboratory  experiments. 
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APPENDIX:  CONDUCTIVE  TEMPERATURE  PROFILE 

Let  the  outer  surfaces  of  the  upper  and  lower  walls  be 
prescribed  at  temperatures  Tl  and  T2,  respectively.  Let  the 
layer  thicknesses  from  below  up  be  dl,  d* ,  and  du.  Use  the 
fluid  layer  thickness  d*  and  the  temperature  difference  A 7'* 
=  Tt-T1  as  the  scales  for  length  and  temperature,  respec¬ 
tively,  so  that  we  note  the  nondimensional  variables  (no  as¬ 
terisk)  T=  flAT*,  z  =  z*/d*,  TX=T */Af*,  T2=f2/AT*,  d, 
=d\  Id  ,  d=  1,  and  du=djd  .  Let  7]  and  74  be  the  nondimen¬ 
sional  temperatures  at  the  inner  upper  and  lower  walls.  The 
constant  dimensional  heat  flux  Q'  in  the  layers  is  given  by 
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(a)  One-Layer,  k=3 


(c)  Three-Layer,  k=3 


(b)  One-Layer,  k=10 


TIME 


FIG.  16.  (Color  online)  Sensor 
temperature  comparison  for  one- 
layer  vs  three-layer,  with  (solid) 
and  without  (dashed)  actuator  lag. 


*  £*Ar*  /  ^ 

Q  —  *  ,  li-l+du  +  dh  (Al) 

d  h 

where  h  is  the  nondimensional  factor  given  in  Eq.  (Al).  The 
temperatures  T3  (at  upper  wall  inner  surface)  and  T4  (at 
lower  wall  inner  surface)  are 


T3  =  T]+dJ(kuh),  T4  =  T2-  djKkjh) .  (A2) 

We  denote  the  thermal  conductivity  ratios  k:=k*/k*  and  ku 
=kjk* .  The  temperature  in  the  three  layers  as  a  function  of  z 
is 


T(z)  = 


T2-(z  +  di)/(kih), 

T4  +  (T3-T4)z, 

J\  +(1  +  du-  z)l(kuh) , 


-d,^z^0, 

0<z«  1, 

1  <  z  =£  1  +du. 

(A3) 


The  effective  Rayleigh  number  is  the  fluid  interface-to- 
interface  Rayleigh  number.  This  number  is  given  by  Raj 
=  a* gAT*(T4-T3)d*3 / v* k* .  We  derive  that 


Raf=l  1-  —  -  —  )Ra, 


kth  k„h 


(A4) 


with  Raj<Ra. 

It  is  helpful  to  show  the  plant  input-output  relationship,  at 
least  for  the  linear  case  in  subcritical  condition.  Consider 
Pr=7.0,  Ra=0.9Rac,  let  the  nominal  wall  thermal  conductiv¬ 


ity  and  diffusivity  values  be  reduced  by  50%  (for  exaggera¬ 
tion).  Consider  a  strong  actuator  lag  with  time  constant  r„ 
=  0.1.  Consider  a  single,  2D  Fourier  mode  of  control  at  wave 
number  of  k  =  3  and  at  A' =10.  The  control  is  the  input  tem¬ 
perature  to  the  plant,  with  amplitude  (Fourier  coefficient)  u 
=  cos(27 r/f)  (here  /=  1).  The  plant  outputs  are  the  tempera¬ 
ture  amplitudes  at  the  nominal  sensor  plane  levels  zJ  =  0.2, 
0.5,  and  0.8.  Figure  16  shows  the  time  responses  of  the  con¬ 
trol.  On  the  left  column,  the  upper  and  lower  panel  corre¬ 
spond  to  the  one-layer  and  three-layer  case,  respectively.  The 
wave  number  of  the  spatial  sinusoidal  disturbance  is  k= 3. 
The  time  responses  in  the  two  panels  are  very  close,  suggest¬ 
ing  that,  at  least  for  the  purely  conductive  case,  the  idealized 
wall  and  finite  wall  boundary  conditions  are  about  the  same. 
The  dashed  curves  (no  lag,  i.e.,  r„ = 0)  are  close  to  the  solid 
curves  (with  lag,  ra=0. 1 )  despite  the  large  actuator  lag  time 
constant.  Note  that  the  peaks  shift  to  the  right  as  zs  is  in¬ 
creased.  However,  for  k=  3  the  temperature  maximum  occurs 
at  the  midplane  rather  than  at  the  sensor  plane  closest  to  the 
actuator.  Similarly,  we  show  the  responses  on  the  right  col¬ 
umn.  The  right  column  corresponds  to  the  larger  wave  num¬ 
ber,  k=  10.  For  the  small-scale  control,  the  temperature  am¬ 
plitude  rapidly  attenuates  when  moved  away  from  the 
actuator  plane.  Unlike  in  the  case  of  k=  3,  the  temperature 
amplitude  drops  monotonically  as  zs  is  increased.  It  suggests 
that  small-scale  control  can  exist  only  in  the  region  close  to 
the  actuator  plane. 
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It  is  shown,  by  direct  numerical  simulations,  that  the  skin-friction  drag  in  a  fully 
developed  channel  can  be  sustained  below  that  corresponding  to  the  laminar  prohle 
when  the  flow  is  subjected  to  surface  blowing  and  suction  in  the  form  of  an  upstream 
travelling  wave.  A  key  mechanism  that  induces  the  sub-laminar  drag  is  the  creation  of 
positive  (negative)  Reynolds  shear  stress  in  the  wall  region,  where  normally  negative 
(positive)  Reynolds  shear  stress  is  expected  given  the  mean  shear.  This  mechanism 
is  contained  in  the  linearized  Navier-Stokes  equations,  thus  allowing  linear  analysis 
of  the  observed  phenomena.  When  applied  to  a  fully  developed  turbulent  channel 
flow,  skin-friction  drag  is  also  significantly  reduced  by  an  upstream  travelling  wave, 
demonstrating  that  the  surface  blowing  and  suction  in  the  form  of  such  a  wave  is 
also  effective  in  fully  developed  turbulent  flows.  Consideration  of  the  energy  budget 
shows  a  possibility  of  net  drag  reduction  in  turbulent  channel  flows  with  the  present 
open-loop  control. 


1.  Introduction 

The  minimum  sustainable  drag  in  a  fully  developed  channel  (pipe)  flow  is  of 
fundamental  importance,  as  it  can  be  used  as  a  basis  for  performance  limitations 
for  various  controllers  designed  to  reduce  skin-friction  drag  in  channel  (pipe)  flow. 
Bewley  (2001,  see  also  Bewley  &  Aamo  2004)  proposed  the  following  conjecture 
(‘Bewley’ s  conjecture’  hereinafter): 

“The  lowest  sustainable  drag  of  an  incompressible  constant  mass-flux  channel  flow,  when  controlled  via 
a  distribution  of  zero-net  mass-flux  blowing/suction  over  the  no-slip  channel  walls,  is  exactly  that  of  the 
laminar  flow." 


This  conjecture  can  be  elucidated  by  a  useful  expression  for  skin-friction  drag  in  fully 
developed  channel  flows  (Fukagata,  Iwamoto  &  Kasagi  2002  and  Bewley  &  Aamo 
2004): 


D-Uf 

d  y 


d U 
d  y 


~2+2Re 


u'v'y  d y. 


(1.1) 


Here,  all  quantities  are  normalized  by  the  centreline  velocity  of  the  laminar  Poiseuille 
flow  (£/c  =  |t/b,  Ub  is  the  bulk  velocity)  and  the  channel  half-height  (5),  U  denotes 
the  mean  velocity,  Re  the  Reynolds  number  based  on  the  laminar  centreline  velocity, 
and  mV  is  the  Reynolds  shear  stress.  In  this  paper,  u,  v,  and  w  denote,  respectively, 
the  velocity  component  in  the  streamwise  ( x ),  wall-normal  (y),  and  spanwise  {z) 
directions,  and  the  prime  denotes  a  fluctuating  quantity.  Note  that  (1.1)  is  valid  for 
all  channel  flows  with  the  same  mass  flux  as  the  base  laminar  flow  (4/3  with  the 
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present  normalization).  The  first  term  on  the  right-hand  side  of  (1.1)  represents  the 
mean  wall-shear  rate  of  laminar  flow  ( U  =  1  — y2),  and  therefore  (1.1)  shows  that 
skin-friction  drag  in  a  channel  flow  consists  of  the  laminar  drag  plus  the  y-weighted 
integral  of  mV.  From  the  viewpoint  of  (1.1),  Bewley’s  conjecture  is  equivalent  to 
saying  that  the  y-weighted  integral  of  u'v',  with  and  without  control  input,  is  always 
positive  in  channel  flows.  For  regular  channel  flows  without  control,  this  is  the 
case,  since  the  Reynolds  shear  stress  in  the  lower  half  of  the  channel  (—  l<y<0) 
is  negative,  while  the  opposite  is  true  in  the  upper  half  of  the  channel  (0  <  y  <  1). 
As  such,  the  skin-friction  drag  in  transitional  and  turbulent  channel  flows  is  higher 
than  that  of  the  corresponding  laminar  flow  with  the  same  mass  flux.  With  a  form 
of  periodic  control,  Bewley  &  Aamo  (2004)  report  that  they  could  not  sustain  the 
Reynolds  shear  stress  necessary  to  yield  drag  below  the  laminar  value.  They  provide 
phenomenological  justification  by  a  Reynolds  analogy  between  convective  momentum 
transport  and  convective  heat  transport,  but  they  left  the  proof  of  the  conjecture  as 
an  open  problem. 

The  objective  of  this  paper  is  to  explore  a  control  input,  in  the  form  of  zero-net- 
mass-flux  blowing  and  suction  on  the  wall,  that  can  sustain  the  Reynolds  shear  stress 
in  such  a  way  that  the  y-weighted  integral  of  u'v'  remains  negative.  In  the  following 
sections,  we  shall  show  that  a  control  input  in  the  form  of  an  upstream  travelling 
wave  indeed  produces  the  Reynolds  shear  stress  that  makes  a  negative  contribution 
to  the  total  drag,  resulting  in  sustained  sub-laminar  drag  in  a  fully  developed  channel 
flow.  It  is  worth  pointing  out  that  although  control  input  in  the  form  of  surface 
blowing  and  suction  used  in  the  present  work  is  easy  to  implement  numerically,  it 
may  not  be  straightforward  to  implement  in  real  flows.  An  alternative  possibility  is 
discussed  in  §  5. 


2.  Linear  analysis 

Recent  studies  have  shown  that  certain  linear  mechanisms  play  important  roles  in 
turbulent  boundary  layers  (for  example,  see  Jimenez  &  Pinelli  1999;  Kim  &  Lim  2000). 
Recognizing  the  role  of  linear  mechanisms,  particularly  that  of  self-sustaining  near¬ 
wall  turbulence  in  turbulent  boundary  layers,  several  investigators  have  successfully 
applied  modern  control  theories  to  develop  optimal  controllers  based  on  the  linearized 
Navier-Stokes  equations  (see,  for  example,  Bewley  2001  and  Kim  2003  and  references 
therein).  The  success  of  these  controllers  when  applied  to  fully  nonlinear  flows  further 
demonstrates  the  usefulness  of  linear  analysis  in  designing  controllers  for  fully  non¬ 
linear  flows.  Here,  we  first  study  the  effect  of  travelling  control  waves  on  the  Reynolds 
shear  stress  by  examining  the  solution  to  the  linearized  Navier-Stokes  equations. 

The  linearized  Navier-Stokes  equation  for  a  two-dimensional  channel  flow  can  be 
written  as 

^.(v=r'(-»t/v=+i„0  +  lv*)o.  in, 


where  v  denotes  the  Fourier-transformed  wall-normal  velocity  perturbation  (i/), 
V2  =  32/3y2  —  a2,  and  a  is  a  wavenumber  in  the  streamwise  direction.  Equation  (2.1) 
can  be  written  in  the  following  state-space  representation  (for  further  details,  see  Kim 
2003): 


dx 

dr 


=  Ax  +  Bu, 


(2.2) 
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Figure  1.  AD  in  steady  state  at  Re  =  2000  and  a  =  0.1:  (a)  downstream  travelling  waves 
(c  >  0);  (b)  upstream  travelling  waves  (c  <  0). 


where  the  vector  x  represents  the  ‘state’  of  the  system  and  consists  of  v  in  Galerkin- 
Chebyshev  space.  The  system  matrix  A  is  related  to  the  linear  operator  in  (2.1),  while 
the  input  matrix  B  and  the  control  u  determine  how  the  control  input  affects  the 
state.  For  any  control  h,  the  solution  to  equation  (2.2)  can  be  found  by: 

x(t)  =  eA'x(0)  +  [  eA{,~T)Bu( r)dr.  (2.3) 

Jo 

We  consider  the  solution  when  control  input  is  introduced  as  surface  blowing  and 
suction  in  the  form  of  a  travelling  wave.  The  initial  objective  of  this  study  was  to 
develop  a  control  strategy  for  viscous  drag  reduction  through  periodic  control  of 
a  turbulent  boundary  layer.  In  the  process  of  optimizing  the  control  input  defined 
at  multiple  wavenumbers,  it  was  observed  that  certain  upstream  travelling  waves 
reduced  the  drag,  while  certain  downstream  travelling  waves  increased  the  drag.  In 
order  to  simplify  our  analysis,  we  consider  travelling  control  waves  consisting  of  a 
single  wavenumber.  Such  a  control  can  be  expressed  in  physical  space  as  boundary 
conditions  for  v : 

uw  =  a  cos(a(x  —  ct)),  (2.4) 

where  a  and  c  denote  the  amplitude  and  wave  speed  of  blowing/suction  on  the  wall, 
respectively.  In  the  present  study,  the  blowing/suction  (2.4)  is  implemented  on  both 
walls  in  varicose  mode,  i.e.  the  upper  and  lower  walls  have  the  blowing/suction  in 
phase  at  the  same  streamwise  locations.  Note  that,  with  stable  systems,  exp(ylf)  — >  0 
as  t— >oo.  Given  control  input  (2.4),  equation  (2.3)  can  be  analytically  solved  for  v  as 
t  — >  oo,  and  the  Reynolds  shear  stress  ( u'v ')  is  obtained  using  the  continuity  equation 
(i  au  +  dv/dy  =  0). 

Figure  1  shows  the  y-weighted  integral  of  u'v'  as  a  function  of  the  wave  speed  for 
Re  =  2000  and  a  =  0.1.  Here,  AH  is  defined  as 


AD 


u'v'  y  dy. 


(2.5) 


Note  that  AD  is  positive  (increased  drag)  with  downstream  travelling  waves  (c  >  0), 
whereas  it  is  negative  (reduced  drag)  with  upstream  travelling  waves  (c  <  0)  except 
near  small  negative  c  (slow  upstream  travelling  wave).  Note  also  that  there  exists  an 
optimal  wave  speed,  where  AD  reaches  its  minimum  (with  negative  c)  and  maximum 
(with  positive  c).  We  mention  in  passing  that  the  large  peak  with  a  positive  c  is 
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related  to  the  most  observable  and  controllable  system  mode.  Also,  although  this 
large  production  of  A D  at  a  certain  c  is  of  no  interest  to  the  present  study,  it  can  be 
useful  in  certain  applications  and  warrants  further  study. 

It  is  apparent  from  the  linear  analysis  that  certain  surface  blowing  and  suction  in 
the  form  of  an  upstream  travelling  wave  can  induce  the  Reynolds  shear  stress  in  such 
a  way  that  the  total  drag  could  be  less  than  that  of  the  laminar  flow.  This  is  contrary 
to  Bewley’ s  conjecture.  Strictly  speaking,  however,  the  linear  analysis  assumes  that 
the  mean  velocity  profile,  and  hence  the  drag,  is  not  affected  by  perturbations  (i.e.  the 
system  matrix  A  is  independent  of  x).  The  real  effect  of  travelling  waves  on  the  drag 
must  be  investigated  by  a  direct  numerical  simulation,  where  the  nonlinear  effects 
are  included.  In  the  following  section,  we  begin  investigating  whether  such  Reynolds 
shear  stress  can  be  sustained  in  nonlinear  flows  with  the  same  control  input. 

Bewley  &  Aamo  (2004)  provided  a  phenomenological  justification  by  an  analogy 
between  convective  momentum  transport  and  convective  heat  transport.  Clearly,  this 
analogy  does  not  hold.  This  failure  of  the  analogy,  however,  is  not  due  to  nonlinearity 
of  momentum  transport  since  the  present  analysis  is  linear.  Rather,  it  is  due  to  dif¬ 
ferences  in  the  linear  equations.  Among  others,  convective  heat  transport  is  governed 
by  a  single  scalar  equation,  whereas  convective  momentum  transport  is  governed  by 
two  equations  coupled  through  the  continuity  equation.  The  present  linear  analysis 
indicates  that  a  certain  flow-field  unsteadiness  induced  by  wall-normal  motion  can 
decelerate  the  momentum  transport  in  the  direction  of  viscous  diffusion,  whereas 
the  heat  transport  is  always  accelerated  by  flow-field  unsteadiness  as  pointed  out  by 
Bewley  &  Aamo  (2004). 


3.  Two-dimensional  channel  flow 

We  will  first  study  the  effect  of  travelling  waves  in  two-dimensional  channel  flows. 
Disturbances  are  finite  but  they  are  two-dimensional,  similar  to  those  considered 
by  Bewley  &  Aamo  (2004).  A  pseudo-spectral  code  similar  to  that  used  by  Kim, 
Moin  &  Moser  (1987)  is  employed.  Simulations  are  conducted  at  Re  =  2000,  and 
the  computational  domain  of  4t 18  x  28  is  used  in  the  streamwise  and  wall-normal 
directions,  respectively,  with  a  32  x  65  grid.  The  mean  pressure  gradient  is  varied  to 
maintain  a  constant  mass  flux,  and  the  total  drag  is  measured  by  averaging  the  mean 
velocity  gradient  on  both  walls. 

Figure  2  shows  the  results  for  a  travelling  wave  of  a  =  0.5  at  c  =  —2.  All  simulations 
start  from  the  laminar  flow  with  no  initial  disturbances,  and  a  steady  state  is  reached 
when  t  >  500  as  shown  in  figure  2(a).  It  is  apparent  that  the  sub-laminar  drag  is 
sustained  for  all  amplitudes  shown.  Figure  2(b)  shows  that  the  results  from  the  linear 
analysis  agree  with  those  from  direct  numerical  simulations  for  small  amplitudes.  As 
expected,  the  nonlinear  results  deviate  from  the  linear  solutions  as  the  amplitude  is 
increased  (figure  2b). 

In  order  to  examine  the  physical  mechanism  responsible  for  drag  reduction  by 
the  upstream  travelling  waves,  we  performed  numerical  simulations  of  ‘channel  flow’ 
without  the  imposed  mean  flow,  i.e.  U  =  0  everywhere  initially  and  the  imposed  mean 
pressure  gradient  was  set  to  zero.  Blowing  and  suction  as  prescribed  in  (2.4)  was 
applied,  and  the  simulations  were  carried  out  until  a  statistically  steady  state  was 
reached.  The  momentum  balance  for  this  flow  is 


d 

dy 


^—u'v'  + 


1  dU\ 
Re  dy  J 


=  0 


or 


—  mV  + 


1  d  U 
Re  dy 


=  C. 


(3.1) 
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Figure  2.  Viscous  drag  in  a  two-dimensional  channel  flow  with  a  =0.5  and  c  =  — 2:  (a)  time 
histories  of  D  for  various  control  amplitudes;  ( b )  AD  in  a  steady  state  as  a  function  of  control 
amplitudes. 
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Figure  3.  Contours  of  v'  (top),  u'  (middle)  and  u'v'  (bottom)  near  the  lower  wall  with 
a  standing  wave.  Contour  levels  for  v',  it'  and  u'v'  are  -0.01-0.01,  -0.02-0.02  and 
-0.0001-0.0001,  respectively  (20  levels). 


The  integration  constant  C  must  be  zero  to  satisfy  the  boundary  conditions.  A 
standing  wave  (c  =  0  in  (2.4))  created  no  net  —u'v'  as  expected  from  the  symmetry 
of  the  problem.  Upstream  travelling  waves  created  positive  u'v'  near  the  lower  wall, 
which  was  balanced  by  a  positive  AU /Ay,  implying  that  U  is  positive.  In  other 
words,  upstream  travelling  waves  induced  a  net  mass  flux  in  the  opposite  direction. 
Downstream  travelling  waves  would  create  the  opposite  effect  due  to  the  apparent 
symmetry.  The  amount  of  induced  mass  flux  was  proportional  to  the  amplitude  of 
the  travelling  wave.  Contours  of  v',  u',  and  u'v'  near  the  lower  wall  with  a  standing 
wave  and  an  upstream-travelling  wave  are  shown  in  figures  3  and  4,  respectively.  It 
is  apparent  that  the  travelling  wave  had  a  different  impact  on  the  phase  of  u'  and 
v' .  Figure  5  illustrates  the  effect  of  standing  and  travelling  waves  on  u'  and  v'  at  a 
location  near  the  lower  wall  (y  =  —0.95)  at  cot  =  2n.  With  a  standing  wave,  the  phase 
of  u',  which  is  related  to  that  of  Av' /Ay  through  the  continuity  equation,  is  90°  out 
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Figure  4.  Contours  of  v'  (top),  u'  (middle)  and  u'v'  (bottom)  near  the  lower  wall  with  an 
upstream  travelling  wave.  Contour  levels  are  the  same  as  those  used  in  figure  3. 
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Figure  5.  Effect  of  standing  and  travelling  waves  on  u!  and  v'  at  y=—  0.95  for  He  =  2000, 

a  =  0.5  and  a  =  0.01 :  (a)  v';  (b)  u'. 


of  phase  as  shown  in  the  figure.  With  an  upstream  travelling  wave,  the  phase  of  v' 
remains  approximately  the  same  as  that  of  the  travelling  wave,  while  that  of  u'  is 
leading  (upstream  direction)  that  of  o'.  This  phase  lead  in  u'  results  in  a  net  positive 
Reynolds  stress  as  can  be  seen  from 

v'  =  ci  cos(ax),  (3.2) 

u  =  cj  sin(ax  +  0),  (3.3) 

u'v '  =  ciC2sin(a.r  +  </>)  cos(ax) 

=  \c\c2  sin  0 .  (3.4) 

With  a  positive  0  as  shown  in  figure  5,  u'v'  is  positive.  With  a  downstream  travelling 
wave,  4>  is  negative,  resulting  in  a  negative  u'v'.  This  phase  shift  in  u'  resulted  in  net 
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Figure  6.  AD  as  a  function  of  the  wave  speed  and  wavenumber  for  a  =0.1. 

positive  (negative)  Reynolds  shear  stress  with  an  upstream  (downstream)  travelling 
wave,  which  in  turn  induced  net  mass  flux  in  the  direction  opposite  to  the  travelling 
wave.  The  effect  of  travelling  waves  in  a  channel  flow  would  be  similar,  although  it 
would  have  been  modified  by  the  presence  of  mean  U  driven  by  the  mean  pressure 
gradient.  Since  an  upstream  travelling  wave  induces  a  net  mass  flux  in  the  streamwise 
direction,  the  amount  of  mass  flux  to  be  driven  by  the  pressure  gradient  in  the  channel 
flow  with  a  fixed  mass  flux  is  reduced,  resulting  in  reduced  drag.  If  we  had  fixed  the 
mean  pressure  gradient  instead,  the  net  mass  flux  in  the  channel  would  have  increased 
due  to  the  positive  mass  flux  induced  by  an  upstream  travelling  wave. 

The  flow  remains  stable  for  large  a  as  the  drag  continues  to  decrease.  However, 
when  the  amplitude  exceeds  a  certain  value,  the  induced  mass  flux  exceeds  the  fixed 
mass  flux.  In  fact,  the  drag  (and  hence  the  pressure  force  required  to  maintain  the 
same  mass  flux)  becomes  negative  as  the  entire  flow  is  driven  by  the  power  required 
to  provide  the  blowing  and  suction,  and  comparison  with  a  channel  flow  is  no  longer 
meaningful.  With  downstream  travelling  waves  (c  >  0),  on  the  other  hand,  the  flow 
becomes  unstable  and  the  drag  increases  drastically  as  a  is  increased. 

Figure  6  shows  drag  variations  as  a  function  of  the  wave  speed  for  different 
streamwise  wavenumbers  with  a  fixed  amplitude  a  =  0.1.  Note  that  AD  in  nonlinear 
simulations  represents  the  change  in  the  total  drag  (see  equation  (1.1)).  The  nonlinear 
results  show  the  same  trend  observed  in  the  linear  solutions:  that  is,  AD  becomes 
larger  for  smaller  a,  and  there  is  an  optimal  wave  speed  that  induces  the  minimum 
AD.  The  optimal  wave  speed  for  nonlinear  simulations  appears  to  be  slightly  more 
negative  (faster  upstream)  than  that  of  the  linear  solutions. 

4.  Turbulent  channel  flow 

In  this  section,  we  investigate  the  effect  of  travelling  waves  in  a  turbulent  channel 
flow.  The  same  code  is  used  to  perform  direct  numerical  simulations  of  a  turbulent 
channel  flow  at  Re  =  2000.  The  computational  domain  is  4nS  x  28  x  4ji<$/3  in  the 
streamwise,  wall-normal  and  spanwise  directions,  respectively,  and  64  x  97  x  64  grid 
points  are  used  in  each  direction.  All  simulations  were  started  from  a  turbulent 
channel  flow  that  had  reached  a  steady  state  without  control  input.  Figure  7  shows 
time  histories  of  D  in  a  turbulent  channel  flow.  An  upstream  travelling  wave  at 
a  =  0.5  with  c  =  —  2  is  applied  on  both  walls  in  varicose  mode,  as  was  the  case  for  the 
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Figure  7.  Time  histories  of  D  in  a  turbulent  channel  flow  for  a  =  0.5  and  c  =  —2. 
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Figure  8.  Distribution  of  total,  viscous,  and  Reynolds  shear  stresses  in  the  upper  half  of  the 
channel  for  different  control  amplitudes.  This  is  for  the  case  of  a  =  0.5  and  c  =  —2. 

two-dimensional  channel  discussed  in  §  3.  The  same  parameters  as  in  two-dimensional 
channel  flow  were  used,  and  they  are  by  no  means  optimal  ones  for  turbulent 
channel  flows.  Similar  to  the  two-dimensional  flows,  a  significant  drag  reduction  is 
obtained,  and  the  reduction  is  greater  for  larger  amplitudes  (about  30  %  and  70  % 
drag  reduction,  respectively,  with  amplitudes  0.1  and  0.25).  Note  that,  with  a  =  0.25, 
the  sustained  drag  is  sub-laminar  even  for  this  three-dimensional  flow.  It  is  worth 
noting,  however,  that  this  amplitude  is  much  larger  than  that  used  in  the  opposition 
control  of  Choi,  Moin  &  Kim  (1994)  and  in  the  LQG  control  of  Lee  et  al.  (2001). 

Figure  8  shows  the  total  shear  stress  for  the  case  of  a  =  0.5  and  c  =  — 2  with 
a  =  0,0.1  and  0.25.  The  total  shear  stress  for  each  case  is  a  straight  line,  indicating 
that  the  flow  has  indeed  reached  a  statistically  steady  state.  Also  shown  in  the  figure 
are  viscous  stress  (dC/dy)  and  Reynolds  shear  stress  (mV).  The  positive  u'v'  in 
uncontrolled  turbulent  flow  (note  that  only  the  upper  half  of  the  channel  is  shown 
in  the  figure,  where  the  Reynolds  shear  stress  is  normally  positive)  is  decreased  by 
the  effect  of  the  upstream  travelling  wave.  For  the  case  of  a  =  0.25,  u'v'  near  the  wall 
becomes  negative,  which  in  turn  results  in  sub-laminar  drag. 
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The  efficiency  of  the  present  control  input  in  the  form  of  a  travelling  wave  is  of 
great  interest.  In  the  present  study,  the  efficiency  p  is  defined  as 

v  -  fdra£  +  PinPut  4 

Po 

Here,  If.  Pdrag  and  PmpM  denote  the  power  required  for  uncontrolled  flow,  the  power 
required  for  drag-reduced  flow,  and  power  input  for  blowing  and  suction,  respectively, 
and  they  can  be  expressed  as 

4  1  4  1  - 

Po  =  2RfD0,  P<ing  =  3  Pinput  =  ( P*  +  Tw)vw-  (4.2) 

where  D0  denotes  the  mean  wall-shear  rate  of  uncontrolled  flow,  and  pw  and  t>w  are 
the  pressure  and  wall-normal  velocity  at  the  wall,  respectively.  Recall  that  the  bulk 
velocity  used  in  defining  the  power  is  2/3  with  the  present  normalization  and  the  total 
drag  is  the  sum  of  the  drag  on  each  wall.  Also  note  that  additional  power  required 
to  account  for  the  friction  loss  associated  with  blowing  and  suction  through  the  wall 
is  not  included  in  (4.1).  For  a  =0.1  and  a  =  0.25,  the  efficiency  was  found  to  be 
0.76  and  0.81,  respectively.  In  other  words,  the  total  power  required  to  have  the  same 
mass  flux  was  only  76%  and  81%  of  the  power  required  in  a  channel  without  control. 
Recall  that  there  was  about  30%  and  70%  drag  reduction  for  a  =  0.1  and  a  =  0.25, 
respectively.  In  comparison,  the  opposition  control  of  Choi  et  al.  (1994)  reduced  drag 
in  a  channel  by  about  30%  and  was  found  to  be  about  0.7.  For  that  closed-loop 
control,  Pinput  was  negligible  compared  to  P0,  and  therefore  the  power  saved  is  directly 
related  to  the  reduced  drag.  For  the  present  open-loop  control,  however,  Pinput  was 
significant,  and  the  power  saved  was  less  (especially  for  high  a)  than  the  saving  due  to 
reduced  drag.  This  was  in  part  due  to  the  large  amplitudes  of  the  travelling  waves.  Its 
low  power  saving  notwithstanding,  it  is  remarkable  that  substantial  drag  reduction 
can  be  achieved  in  a  turbulent  channel  flow  with  such  a  simple  open-loop  control. 

5.  Concluding  remarks 

Motivated  by  Bewley’s  conjecture  (Bewley  2001),  we  investigated  the  possibility  of 
achieving  sub-laminar  drag  in  a  fully  developed  channel  flow.  Sustainable  sub-laminar 
drag  was  obtained  when  the  flow  was  subjected  to  blowing  and  suction  at  the  wall  in 
the  form  of  an  upstream  travelling  wave.  It  was  found,  both  from  linear  analysis  and 
nonlinear  simulations,  that  certain  upstream  travelling  waves  induce  the  Reynolds 
shear  stress  in  such  a  way  that  it  makes  a  negative  contribution  to  the  total  viscous 
drag.  This  was  the  case  not  only  for  the  two-dimensional  channel  flow  considered 
by  Bewley  &  Aamo  (2004),  but  also  for  a  three-dimensional  turbulent  channel  flow. 
Note  that  the  Reynolds  shear  stress  induced  by  upstream  travelling  waves  reduces  the 
production  of  kinetic  energy,  and  therefore  the  flow  remains  stable  for  large-amplitude 
upstream  travelling  waves. 

Downstream  travelling  waves,  on  the  other  hand,  increase  the  drag.  The  linear 
analysis  shows  that  at  certain  wave  speeds  the  drag  increase  is  dramatic.  There  are 
certain  applications  where  the  increase  in  Reynolds  shear  stress  could  be  desirable. 
For  example,  an  optimized  downstream  travelling  wave  could  be  used  to  prevent  or 
delay  separation  in  turbulent  boundary  layers  subject  to  a  strong  pressure  gradient 
(e.g.  in  a  diffuser  or  flow  over  an  airfoil).  The  optimal  use  of  downstream  travelling 
waves  is  something  that  warrants  further  study. 

The  open-loop  control  presented  was  discovered  serendipitously  while  exploring 
periodic  control  of  turbulent  boundary  layers.  Unlike  the  cyclic  application  of  a 
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pressure  feedback  control  (turning  on  and  off  their  ‘win-win’  mechanism)  by  Bewley  & 
Aamo  (2004),  our  periodic  optimization  specifically  tasked  the  wall-normal  surface 
blowing/suction  control  to  return  the  flow  state  back  to  its  initial  condition.  Both 
the  time  period  of  control  and  initial  state  were  part  of  the  optimization  (see  Speyer 
1996).  Numerically  calculated  gradients  were  used  to  find  a  control  history,  an  initial 
state,  and  a  time  period  that  minimized  viscous  drag.  The  two-dimensional  travelling 
wave  was  discovered  as  we  simplified  the  optimization  to  make  the  problem  more 
tenable.  The  control  presented  in  this  paper  is  not  an  optimal  solution;  our  purpose 
was  to  investigate  whether  a  sub-laminar  drag  can  be  sustained  in  a  fully  developed 
channel  flow.  In  the  light  of  this  new  finding,  optimization  is  underway.  In  particular, 
we  plan  to  explore  control  input  in  the  form  of  spanwise  and  obliquely  travelling 
waves  as  well  as  constructing  periodic  regulators  to  form  closed-loop  solutions. 

Finally,  the  current  control  scheme,  consisting  of  surface  blowing  and  suction  in  the 
form  of  travelling  waves,  is  mathematically  simple  (and  hence  easy  to  implement  in 
numerical  experiments),  yet  it  may  not  be  straightforward  to  implement  in  real  flows. 
For  example,  additional  space  is  required  to  retain  the  flow  through  the  walls,  and 
the  additional  hardware  required  to  control  blowing  and  suction  may  be  complicated. 
However,  a  moving  surface  with  wavy  motion  would  produce  a  similar  effect,  since 
wavy  walls  with  small  amplitudes  can  be  approximated  by  surface  blowing  and 
suction.  We  plan  to  perform  simulations  over  moving  wavy  walls. 
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